Mid-Cayman Rise biomass analysis
Sarah Hu
2/3/2021
RMarkdown history: 1. Feb 2021 started github repo (https://github.com/shu251/midcayman-rise-microeuk) 2. continue to commit changes and as you include knitr and commit changes, those will be pushed to the git repo 3. eventually set up the gitpages under settings 4 https://resources.github.com/whitepapers/github-and-rstudio/
1 Mid-Cayman Rise grazing
Summary of cruise operations
Below is the compilation of all microscopy data
1.1 Set up working environment
GitHub repo hosts all raw data to reproduce analysis and results. See input-data or clone the repo.
library(tidyverse); library(cowplot); library(broom)
getwd()## [1] "/Users/sarahhu/Desktop/Projects/Mid-Cayman_Rise/midcayman-rise-microeuk"2 Protistan cell concentration
Import eukaryotic cell count data from grazing experiments. In this section, we will calculate cells per ml from raw counts (Field of view, etc.) and use to estimate protist cell concentration. These will be used below in grazing experiment calculations.
counts <- read.delim("input-data/euk-counts-compiled.txt",
blank.lines.skip = FALSE,
na.strings = c("", "NA"),
stringsAsFactors = FALSE) # Import
counts[is.na(counts)] <- 0 # Change blanks to zeroesRaw data table collected during microscopy count process. Below code reviews the structure of this raw data and updates column headers to be more ‘R’ friendly.
colnames(counts) <-c("DATE", "SAMPLE", "EXPID", "VOL", "MAG", "FOV", "nanoNoFLP", "microNoFLP", "nanoFLP", "microFLP", "NOTES", "DateCompiled"); colnames(counts)## [1] "DATE" "SAMPLE" "EXPID" "VOL" "MAG"
## [6] "FOV" "nanoNoFLP" "microNoFLP" "nanoFLP" "microFLP"
## [11] "NOTES" "DateCompiled"# Column 9 and 10 list the total number of FLP inside each, so # of occurences
# Notes if count was not countable
# unique(counts$NOTES)To count occurence and number of FLP ingested by eukaryotic cells, the number of FLPs ingested was tallied and comma separated for multiple eukaryotic cells with FLP. These values need to separated and counted as 1 eukaryotic cell each, but retain the number of FLP per cell.
counts_occur <- counts %>%
filter(NOTES != "Not countable") %>%
mutate(nanoFLP_occur = as.numeric(str_count(nanoFLP, "[1-9]\\d*")), #Count number of euk cells observed with FLPs
microFLP_occur = as.numeric(str_count(microFLP, "[1-9]\\d*")),
nanoTOTAL = as.numeric(nanoNoFLP) + nanoFLP_occur, #Add number of euk cells with FLPs to those without for total number of euk cells
microTOTAL = as.numeric(microNoFLP) + microFLP_occur,
euksTOTAL = nanoTOTAL + microTOTAL) %>%
data.frame2.1 Microscopy cells per ml (euk)
Input data are the raw microscopy counts by FOV. Code below calculations cells/ml based on these values. Additionaly, variance and standard deviation are also calculated. Eukaryotic cells were also classified by size, where micro equates to >20um and nano is <20um. All counts were done at 100x magnification, confirm this:
unique(counts_occur$MAG)## [1] 100# head(counts_occur)Calculate cell concentration (cells/ml).
counts_cellsml_all <- counts_occur %>%
group_by(SAMPLE, EXPID, VOL) %>% #Calculate averages by sample
summarise(totalFOV = n(), # Count total FOV counted
nanoAvg = sum(nanoTOTAL)/totalFOV, #Average per FOV
nanoVar = var(nanoTOTAL), #Variance
nanoSd = (2*(sqrt(nanoVar))), #Standard deviation
microAvg = sum(microTOTAL)/totalFOV, ## Repeat for microeuks
microVar = var(microTOTAL),
microSd = (2*(sqrt(microVar))),
euksAvg = sum(euksTOTAL)/totalFOV, ## Repeat for total cell count
euksVar = var(euksTOTAL),
euksSd = (2*(sqrt(euksVar))),
.groups = 'drop_last') %>%
# Calculate cells/ml based on magnification (at x100, 0.01 is vol of grid), volume filtered (VOL), dilution factor (0.9), and area of counting grid (for Huber lab scope, it is 283.385):
mutate(nanoCONC = ((nanoAvg * 283.385)/(VOL * 0.01 * 0.9)),
microCONC = ((microAvg * 283.385)/(VOL * 0.01 * 0.9)),
eukCONC = ((euksAvg * 283.385)/(VOL * 0.01 * 0.9))
) %>%
# left_join(expmeta) %>%
separate(SAMPLE, c("Site", "Name"), sep = "-", remove = FALSE) %>%
separate(EXPID, c("TimePoint", "Replicate"), sep = "-", remove = FALSE) %>%
data.framehead(counts_cellsml_all)## SAMPLE Site Name EXPID TimePoint Replicate VOL
## 1 Piccard-LotsOShrimp Piccard LotsOShrimp T0-Rep3 T0 Rep3 50
## 2 Piccard-LotsOShrimp Piccard LotsOShrimp T15-Rep3 T15 Rep3 50
## 3 Piccard-LotsOShrimp Piccard LotsOShrimp T20-Rep3 T20 Rep3 50
## 4 Piccard-LotsOShrimp Piccard LotsOShrimp T40-Rep3 T40 Rep3 40
## 5 Piccard-Plume Piccard Plume T0-Rep1 T0 Rep1 150
## 6 Piccard-Plume Piccard Plume T0-Rep2 T0 Rep2 150
## totalFOV nanoAvg nanoVar nanoSd microAvg microVar microSd
## 1 30 0.3666667 0.3781609 1.2298958 0.00000000 0.00000000 0.0000000
## 2 30 0.3000000 0.2172414 0.9321832 0.06666667 0.06436782 0.5074163
## 3 30 0.3333333 0.2988506 1.0933445 0.00000000 0.00000000 0.0000000
## 4 30 0.1666667 0.1436782 0.7580980 0.00000000 0.00000000 0.0000000
## 5 30 0.3000000 0.2862069 1.0699662 0.03333333 0.03333333 0.3651484
## 6 30 0.2000000 0.1655172 0.8136762 0.06666667 0.06436782 0.5074163
## euksAvg euksVar euksSd nanoCONC microCONC eukCONC
## 1 0.3666667 0.3781609 1.229896 230.90630 0.00000 230.90630
## 2 0.3666667 0.2402299 0.980265 188.92333 41.98296 230.90630
## 3 0.3333333 0.2988506 1.093345 209.91481 0.00000 209.91481
## 4 0.1666667 0.1436782 0.758098 131.19676 0.00000 131.19676
## 5 0.3333333 0.3678161 1.212957 62.97444 6.99716 69.97160
## 6 0.2666667 0.2712644 1.041661 41.98296 13.99432 55.97728Replicates belong to the same experiment for either Bag or IGT incubation. Below, modify these names and label new column with bag or igt. And create an average across replicates.
Average cells/ml across replicates, pivot to long format
counts_cellsml_avg <- counts_cellsml_all %>%
select(Site, Name, TimePoint, Replicate, nanoCONC, microCONC, eukCONC) %>%
mutate(EXP_TYPE = case_when(
grepl("IGT", Replicate) ~ "IGT",
grepl("Rep", Replicate) ~ "Bag"
)) %>%
mutate(IGT_REP = case_when(
EXP_TYPE == "IGT" ~ Replicate,
EXP_TYPE == "Bag" ~ "Bag")) %>%
select(-Replicate) %>%
pivot_longer(cols = ends_with("CONC"), names_to = "VARIABLE", values_to = "CONCENTRATION") %>%
group_by(Site, Name, TimePoint, EXP_TYPE, IGT_REP, VARIABLE) %>%
# Calculate mean, variance, SD, min, and max
summarise(MEAN = mean(CONCENTRATION),
VAR = var(CONCENTRATION),
SD = sd(CONCENTRATION),
SEM =(sd(CONCENTRATION)/sqrt(length(CONCENTRATION))),
MIN = min(CONCENTRATION),
MAX = max(CONCENTRATION),
.groups = 'drop_last') %>%
data.frame
head(counts_cellsml_avg)## Site Name TimePoint EXP_TYPE IGT_REP VARIABLE MEAN VAR SD SEM
## 1 Piccard LotsOShrimp T0 Bag Bag eukCONC 230.90630 NA NA NA
## 2 Piccard LotsOShrimp T0 Bag Bag microCONC 0.00000 NA NA NA
## 3 Piccard LotsOShrimp T0 Bag Bag nanoCONC 230.90630 NA NA NA
## 4 Piccard LotsOShrimp T15 Bag Bag eukCONC 230.90630 NA NA NA
## 5 Piccard LotsOShrimp T15 Bag Bag microCONC 41.98296 NA NA NA
## 6 Piccard LotsOShrimp T15 Bag Bag nanoCONC 188.92333 NA NA NA
## MIN MAX
## 1 230.90630 230.90630
## 2 0.00000 0.00000
## 3 230.90630 230.90630
## 4 230.90630 230.90630
## 5 41.98296 41.98296
## 6 188.92333 188.92333# View(counts_cellsml_avg)
# NOTES on calculations:
# VAR = takes the sum of the squares of each value's deviation from the mean and divides by the number of such values minus one. This differs from the calculation of variance across an entire population in that the latter divides by the size of the dataset without subtracting one.
# SD = standard deviation of all values
# SEM = standard deviation of sampling distribution; standard deviation divided by the square root of the sample size.# Euk cell counts, all and averaged
head(counts_cellsml_all)## SAMPLE Site Name EXPID TimePoint Replicate VOL
## 1 Piccard-LotsOShrimp Piccard LotsOShrimp T0-Rep3 T0 Rep3 50
## 2 Piccard-LotsOShrimp Piccard LotsOShrimp T15-Rep3 T15 Rep3 50
## 3 Piccard-LotsOShrimp Piccard LotsOShrimp T20-Rep3 T20 Rep3 50
## 4 Piccard-LotsOShrimp Piccard LotsOShrimp T40-Rep3 T40 Rep3 40
## 5 Piccard-Plume Piccard Plume T0-Rep1 T0 Rep1 150
## 6 Piccard-Plume Piccard Plume T0-Rep2 T0 Rep2 150
## totalFOV nanoAvg nanoVar nanoSd microAvg microVar microSd
## 1 30 0.3666667 0.3781609 1.2298958 0.00000000 0.00000000 0.0000000
## 2 30 0.3000000 0.2172414 0.9321832 0.06666667 0.06436782 0.5074163
## 3 30 0.3333333 0.2988506 1.0933445 0.00000000 0.00000000 0.0000000
## 4 30 0.1666667 0.1436782 0.7580980 0.00000000 0.00000000 0.0000000
## 5 30 0.3000000 0.2862069 1.0699662 0.03333333 0.03333333 0.3651484
## 6 30 0.2000000 0.1655172 0.8136762 0.06666667 0.06436782 0.5074163
## euksAvg euksVar euksSd nanoCONC microCONC eukCONC
## 1 0.3666667 0.3781609 1.229896 230.90630 0.00000 230.90630
## 2 0.3666667 0.2402299 0.980265 188.92333 41.98296 230.90630
## 3 0.3333333 0.2988506 1.093345 209.91481 0.00000 209.91481
## 4 0.1666667 0.1436782 0.758098 131.19676 0.00000 131.19676
## 5 0.3333333 0.3678161 1.212957 62.97444 6.99716 69.97160
## 6 0.2666667 0.2712644 1.041661 41.98296 13.99432 55.97728head(counts_cellsml_avg)## Site Name TimePoint EXP_TYPE IGT_REP VARIABLE MEAN VAR SD SEM
## 1 Piccard LotsOShrimp T0 Bag Bag eukCONC 230.90630 NA NA NA
## 2 Piccard LotsOShrimp T0 Bag Bag microCONC 0.00000 NA NA NA
## 3 Piccard LotsOShrimp T0 Bag Bag nanoCONC 230.90630 NA NA NA
## 4 Piccard LotsOShrimp T15 Bag Bag eukCONC 230.90630 NA NA NA
## 5 Piccard LotsOShrimp T15 Bag Bag microCONC 41.98296 NA NA NA
## 6 Piccard LotsOShrimp T15 Bag Bag nanoCONC 188.92333 NA NA NA
## MIN MAX
## 1 230.90630 230.90630
## 2 0.00000 0.00000
## 3 230.90630 230.90630
## 4 230.90630 230.90630
## 5 41.98296 41.98296
## 6 188.92333 188.92333Parse experiment type information
# Convert to long format and add column that reports IGT vs bag experiment
plot_euk_conc <- counts_cellsml_all %>%
select(Site, Name, TimePoint, Replicate, ends_with("CONC")) %>%
mutate(EXP_TYPE = case_when(
grepl("IGT", Replicate) ~ "IGT",
grepl("Rep", Replicate) ~ "Bag"
)) %>%
pivot_longer(cols = ends_with("CONC"), names_to = "VARIABLE", values_to = "CONCENTRATION") %>%
data.frame# head(plot_euk_conc)
unique(plot_euk_conc$Name)## [1] "LotsOShrimp" "Plume" "Shrimpocalypse" "BSW"
## [5] "MustardStand" "OMT" "Rav2" "ShrimpHole"
## [9] "X18"vent_ids <- c("BSW","Plume", "Shrimpocalypse", "LotsOShrimp", "X18", "OMT", "Rav2", "MustardStand", "ShrimpHole")
vent_fullname <- c("Background","Plume", "Shrimpocalypse", "Lots 'O Shrimp", "X-18", "Old Man Tree", "Ravelin #2", "Mustard Stand", "Shrimp Hole")
site_ids <- c("VD", "Piccard")
site_fullname <- c("Von Damm", "Piccard")
plot_euk_conc$SiteOrder <- factor(plot_euk_conc$Site, levels = site_ids, labels = site_fullname)
plot_euk_conc$NameOrder <- factor(plot_euk_conc$Name, levels = vent_ids, labels = vent_fullname)Box plot to report all eukaryote cell counts.
# Code for plot
conc_boxplot <- ggplot(plot_euk_conc, aes(x = NameOrder,
y = CONCENTRATION,
group = NameOrder,
fill = VARIABLE,
shape = EXP_TYPE)) +
geom_boxplot() +
# Do not color by time point
geom_jitter(color = "black", size = 2, aes(fill = VARIABLE,
shape = EXP_TYPE)) +
scale_shape_manual(values = c(21,24)) +
scale_fill_manual(values = c("#e7298a", "#fcbba1", "#c6dbef")) +
coord_flip() +
scale_y_log10() +
# scale_y_log10(limits = c(10,1000), expand = c(0, 0)) +
facet_grid(SiteOrder ~ EXP_TYPE, space = "free", scale = "free") +
theme_bw() +
theme(axis.text.x = element_text(angle = 0, h = 1, vjust = 1),
strip.background = element_blank(),
legend.position = "right",
legend.title = element_blank()) +
labs(x = "", y = bquote("Eukaryote cells "~mL^-1),
title = "Distribution of all eukaryotic cell counts")Eukaryote cell concentration (cells/ml) are lower in the background and plume samples compared to vent sites. ~300 cells/ml in background and plume compared to ~1000 cells per ml at the vent sites. These values are also consistent between each vent site (Von Damm and Piccard) and between Bag and IGT samples.
Boxplot represents the median (line in box) and the 1st and 3rd quartiles in the lower and upper hinges, respectively (25th and 75th percentiles). Black data points are outliers from the boxplot. Upper and lower whiskers represent the 1.5 * interquartile ranges. Pink data points are the values contributing to the boxplot (individial counts across replicates and time points.)
conc_boxplot## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Transformation introduced infinite values in continuous y-axis## Warning: Removed 39 rows containing non-finite values (stat_boxplot).
eukCONC is the sum of micro and nano. Because there was a discrepency between the micro and nano cell counts, we plan to combine for most of the analysis. Here we show that the cell concentration across replicate samples was similar throughout experiments. And that the bag versus IGT experiment results were within range of one another.
Subset to plot from T0 only.
vent_ids <- c("BSW","Plume", "Shrimpocalypse", "LotsOShrimp", "X18", "OMT", "Rav2", "MustardStand", "ShrimpHole")
vent_fullname <- c("Background","Plume", "Shrimpocalypse", "Lots 'O Shrimp", "X-18", "Old Man Tree", "Ravelin #2", "Mustard Stand", "Shrimp Hole")
plot_euk_format <- plot_euk_conc %>%
filter(TimePoint == "T0" & (VARIABLE == "eukCONC")) %>%
group_by(SiteOrder, NameOrder, TimePoint, EXP_TYPE, VARIABLE) %>%
summarise(avg_conc = mean(CONCENTRATION),
SEM_conc = (sd(CONCENTRATION)/sqrt(length(CONCENTRATION))),
.groups = "rowwise") %>%
unite(EXPERIMENT, SiteOrder, NameOrder, EXP_TYPE, remove = FALSE) %>%
data.frame
# Factor
plot_euk_format$Site_Order <- factor(plot_euk_format$SiteOrder, levels = site_fullname, labels = site_fullname)
# View(plot_euk_format)
euk_plot <- ggplot(plot_euk_format, aes(x = NameOrder, y = avg_conc, fill = Site_Order)) +
geom_errorbar(aes(ymax = (avg_conc + SEM_conc), ymin = (avg_conc - SEM_conc)), width = 0.2) +
geom_point(aes(fill = Site_Order), color = "black", stat = "identity", size = 3, shape = 23) +
facet_grid(.~ Site_Order, space = "free", scales = "free") +
scale_fill_manual(values = c("#1c9099", "#de2d26")) +
theme_minimal() +
theme(panel.grid.major = element_line(), panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
axis.text.x = element_text(color="black", size = 12,
angle = 45, hjust = 1, vjust = 1),
axis.text.y = element_text(color="black", size = 12),
axis.title =element_text(color="black", size = 12),
axis.ticks = element_line(),
strip.text =element_blank(), legend.title = element_blank()) +
labs(x = "", y = bquote("Eukaryote cells "~mL^-1),
title = "")euk_plot
3 Prokaryote cell concentration
Import counts from DAPI counts of bacteria and archaea.
prok <- read.delim("input-data/prokINSITU-counts-compiled.txt")
insitu_proks <- prok %>%
separate(SAMPLE, c("Site", "Name"), sep = "-", remove = FALSE) %>%
group_by(SAMPLE, Site, Name) %>%
summarise(MEAN = mean(CELLML),
SD = sd(CELLML),
SEM = (sd(CELLML)/sqrt(length(CELLML))),
.groups = "rowwise") %>%
data.frame
# head(insitu_proks)3.1 Microscopy cells per ml (prok)
Plot in situ prokaryote values. Average across dives where applicable.
insitu_proks$Name_order <- factor(insitu_proks$Name, levels = c("BSW", "Plume", "Quakeplume", "Shrimpocalypse", "LotsOShrimp", "X18", "OMT", "Rav2", "MustardStand", "ShrimpHole", "HotChimlet1", "ShrimpGulley", "SouthofHotChimlet", "SouthofLungSnack", "ArrowLoop", "Bartizan", "Rav1"), labels = c("Background","Plume", "Quakeplume", "Shrimpocalypse", "Lots 'O Shrimp", "X-18", "Old Man Tree", "Ravelin #2", "Mustard Stand", "Shrimp Hole", "Hot Chimlet #1", "Shrimp Gulley", "South of Hot Chimlet", "South of LungSnack", "Arrow Loop", "Bartizan", "Ravelin #1"))
site_ids <- c("VD", "Piccard")
site_fullname <- c("Von Damm", "Piccard")
insitu_proks$Site_order <- factor(insitu_proks$Site, levels = site_ids, labels = site_fullname)
# unique(insitu_proks$Name)
prok_plot <- ggplot(insitu_proks, aes(x = Name_order, y = MEAN)) +
geom_errorbar(aes(ymax = (MEAN + SEM), ymin = (MEAN - SEM)), width = 0.2) +
geom_point(stat = "identity", shape = 23, aes(fill = Site), size = 3) +
facet_grid(.~ Site_order, space = "free", scales = "free") +
scale_fill_manual(values = c("#de2d26", "#1c9099")) +
labs(y = bquote("Prokaryote cells "~mL^-1), x = "", title = "") +
scale_y_log10() +
theme_minimal() +
theme(panel.grid.major = element_line(), panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
axis.text.x = element_text(color="black", size = 12,
angle = 45, hjust = 1, vjust = 1),
axis.text.y = element_text(color="black", size = 12),
axis.title =element_text(color="black", size = 12),
axis.ticks = element_line(),
strip.text =element_blank(), legend.title = element_blank())prok_plot
3.2 Comparison with previous MCR
Can include years 2012 and 2013, seeking clarification on these data. ## Prokaryote biomass-literature
Work on this! get FS xx vlaues from Julie
4 Determine grazing effect
Calculate FLP per eukaryotic cell over time. Goal is to make these calculations and then determine best fit line. Slope of best fit line is the grazing rate. Need to take into account euk cells with FLPs and then the euk cells withOUT FLPs, these will be zeroes to take into account for FLPs/euk averages.
Retain previously generated dataframes:
# Euk cell counts, all and averaged
head(counts_cellsml_all)## SAMPLE Site Name EXPID TimePoint Replicate VOL
## 1 Piccard-LotsOShrimp Piccard LotsOShrimp T0-Rep3 T0 Rep3 50
## 2 Piccard-LotsOShrimp Piccard LotsOShrimp T15-Rep3 T15 Rep3 50
## 3 Piccard-LotsOShrimp Piccard LotsOShrimp T20-Rep3 T20 Rep3 50
## 4 Piccard-LotsOShrimp Piccard LotsOShrimp T40-Rep3 T40 Rep3 40
## 5 Piccard-Plume Piccard Plume T0-Rep1 T0 Rep1 150
## 6 Piccard-Plume Piccard Plume T0-Rep2 T0 Rep2 150
## totalFOV nanoAvg nanoVar nanoSd microAvg microVar microSd
## 1 30 0.3666667 0.3781609 1.2298958 0.00000000 0.00000000 0.0000000
## 2 30 0.3000000 0.2172414 0.9321832 0.06666667 0.06436782 0.5074163
## 3 30 0.3333333 0.2988506 1.0933445 0.00000000 0.00000000 0.0000000
## 4 30 0.1666667 0.1436782 0.7580980 0.00000000 0.00000000 0.0000000
## 5 30 0.3000000 0.2862069 1.0699662 0.03333333 0.03333333 0.3651484
## 6 30 0.2000000 0.1655172 0.8136762 0.06666667 0.06436782 0.5074163
## euksAvg euksVar euksSd nanoCONC microCONC eukCONC
## 1 0.3666667 0.3781609 1.229896 230.90630 0.00000 230.90630
## 2 0.3666667 0.2402299 0.980265 188.92333 41.98296 230.90630
## 3 0.3333333 0.2988506 1.093345 209.91481 0.00000 209.91481
## 4 0.1666667 0.1436782 0.758098 131.19676 0.00000 131.19676
## 5 0.3333333 0.3678161 1.212957 62.97444 6.99716 69.97160
## 6 0.2666667 0.2712644 1.041661 41.98296 13.99432 55.97728head(counts_cellsml_avg)## Site Name TimePoint EXP_TYPE IGT_REP VARIABLE MEAN VAR SD SEM
## 1 Piccard LotsOShrimp T0 Bag Bag eukCONC 230.90630 NA NA NA
## 2 Piccard LotsOShrimp T0 Bag Bag microCONC 0.00000 NA NA NA
## 3 Piccard LotsOShrimp T0 Bag Bag nanoCONC 230.90630 NA NA NA
## 4 Piccard LotsOShrimp T15 Bag Bag eukCONC 230.90630 NA NA NA
## 5 Piccard LotsOShrimp T15 Bag Bag microCONC 41.98296 NA NA NA
## 6 Piccard LotsOShrimp T15 Bag Bag nanoCONC 188.92333 NA NA NA
## MIN MAX
## 1 230.90630 230.90630
## 2 0.00000 0.00000
## 3 230.90630 230.90630
## 4 230.90630 230.90630
## 5 41.98296 41.98296
## 6 188.92333 188.92333head(counts_occur)## DATE SAMPLE EXPID VOL MAG FOV nanoNoFLP microNoFLP nanoFLP microFLP
## 1 4/2/20 VD-Plume T0-Rep1 150 100 1 0 0 0 0
## 2 4/2/20 VD-Plume T0-Rep1 150 100 2 0 0 0 0
## 3 4/2/20 VD-Plume T0-Rep1 150 100 3 3 0 0 0
## 4 4/2/20 VD-Plume T0-Rep1 150 100 4 1 0 4 0
## 5 4/2/20 VD-Plume T0-Rep1 150 100 5 0 1 0 0
## 6 4/2/20 VD-Plume T0-Rep1 150 100 6 0 0 0 0
## NOTES DateCompiled nanoFLP_occur microFLP_occur nanoTOTAL
## 1 0 6/7/20 0 0 0
## 2 0 6/7/20 0 0 0
## 3 0 6/7/20 0 0 3
## 4 0 6/7/20 1 0 2
## 5 elongated cell? 6/7/20 0 0 0
## 6 0 6/7/20 0 0 0
## microTOTAL euksTOTAL
## 1 0 0
## 2 0 0
## 3 0 3
## 4 0 2
## 5 1 1
## 6 0 0Supplementary plot showing the concentration of eukaryotic cells over time for each experiment.
# head(counts_cellsml_avg)
vent_ids <- c("BSW","Plume", "Shrimpocalypse", "LotsOShrimp", "X18", "OMT", "Rav2", "MustardStand", "ShrimpHole")
vent_fullname <- c("Background","Plume", "Shrimpocalypse", "Lots 'O Shrimp", "X-18", "Old Man Tree", "Ravelin #2", "Mustard Stand", "Shrimp Hole")
site_ids <- c("VD", "Piccard")
site_fullname <- c("Von Damm", "Piccard")
counts_cellsml_avg$SiteOrder <- factor(counts_cellsml_avg$Site, levels = site_ids, labels = site_fullname)
counts_cellsml_avg$NameOrder <- factor(counts_cellsml_avg$Name, levels = vent_ids, labels = vent_fullname)
# Plot trend line of euk cell count for all experiments
counts_cellsml_avg %>%
filter(VARIABLE == "eukCONC") %>%
unite("Experiment", NameOrder, IGT_REP, sep = "-", remove = FALSE) %>%
ggplot(aes(x = TimePoint, y = MEAN, shape = EXP_TYPE, fill = NameOrder)) +
geom_path(aes(group = Experiment)) +
# geom_errorbar(aes(ymax = (MEAN + SD), ymin = (MEAN - SD)), width = 0.2) +
geom_errorbar(aes(ymax = (MEAN + SEM), ymin = (MEAN - SEM)), width = 0.2) +
geom_point(stat = "identity", size = 2, aes(shape = EXP_TYPE)) +
scale_shape_manual(values = c(21, 24)) +
scale_fill_brewer(palette = "Paired") +
scale_y_log10() +
facet_wrap(SiteOrder ~ EXP_TYPE, scales = "free") +
theme_classic() + theme(strip.background = element_blank(),
legend.title = element_blank(),
title = element_text(size = 7, face = "bold"),
axis.title = element_text(size = 9)) +
labs(title = "Total euk cell counts for each experiment", y = bquote("Average eukaryote cells "~mL^-1), x = "Time point") +
guides(fill=guide_legend(override.aes=list(shape=21)))
## Determine FLP per euk Isolate the euk cell counts with FLPs, these are comma separated and must be separated into rows. Using counts_occur data frame from above.
# Select nano and micro counts with FLPs
counts_sepflp <- counts_occur %>%
filter(!NOTES == "Discard") %>%
filter(!(NOTES == "DTAF stain prevented counts of FLP, Euks only")) %>%
select(DATE, SAMPLE, EXPID, VOL, MAG, FOV, nanoFLP, microFLP) %>%
separate_rows(microFLP, sep = ",", convert = TRUE) %>%
separate_rows(nanoFLP, sep = ",", convert = TRUE) %>%
replace_na(list(microFLP = 0, nanoFLP = 0)) %>% #Replace NAs with zero
data.frame
## Check, see FOV 23, separated into rows.
# View(counts_sepflp %>%
# filter(SAMPLE == "VD-Rav2" & EXPID == "T10-Rep1"))
# View(counts_occur %>%
# filter(SAMPLE == "VD-Rav2" & EXPID == "T10-Rep1"))Subset counts that are greater than 0, so only euk cells observed to have FLPs. Create new column that calculates FLP per Euk - by dividing by 1.
counts_flp <- counts_sepflp %>%
select(SAMPLE, EXPID, nano_size = nanoFLP, micro_size = microFLP) %>%
pivot_longer(cols = ends_with("_size"), names_to = "SizeFrac", values_to = "num_of_FLP") %>%
filter(num_of_FLP > 0) %>%
separate(SAMPLE, c("Site", "Name"), sep = "-", remove = FALSE) %>%
separate(EXPID, c("TimePoint", "Replicate"), sep = "-", remove = FALSE) %>%
mutate(EXP_TYPE = case_when(
grepl("IGT", Replicate) ~ "IGT",
grepl("Rep", Replicate) ~ "Bag"
)) %>%
mutate(IGT_REP = case_when(
EXP_TYPE == "IGT" ~ Replicate,
EXP_TYPE == "Bag" ~ "Bag")) %>%
group_by(SAMPLE, EXPID, EXP_TYPE, IGT_REP, SizeFrac) %>%
summarise(total_FLP = sum(num_of_FLP),
total_euks_wflp = n(),
.groups = "rowwise") %>%
data.frameOUTPUT COLUMNS: (1) total_FLP = sum of FLPs found inside a euk cell (2) total_euks_wflp = number of euks counted with ingested FLP
Repeat above operation for euk cells without any FLP. Here, subset total number of observations where there was a euk cell without FLP. These need to be counted as euk cell without an FLP. > Below code repeats process and compiles with other FLP/euk cell data.
counts_flp_compiled <- counts_occur %>%
filter(!(NOTES == "Discard")) %>% #Discard bad counts
filter(!(NOTES == "DTAF stain prevented counts of FLP, Euks only")) %>%
type.convert(as.is = TRUE) %>% #modify str() for columns
select(SAMPLE, EXPID, nano_size = nanoNoFLP, micro_size = microNoFLP) %>% #select non flp
pivot_longer(cols = ends_with("_size"), names_to = "SizeFrac", values_to = "num_of_euks") %>%
separate(SAMPLE, c("Site", "Name"), sep = "-", remove = FALSE) %>%
separate(EXPID, c("TimePoint", "Replicate"), sep = "-", remove = FALSE) %>%
mutate(EXP_TYPE = case_when(
grepl("IGT", Replicate) ~ "IGT",
grepl("Rep", Replicate) ~ "Bag"
)) %>%
mutate(IGT_REP = case_when(
EXP_TYPE == "IGT" ~ Replicate,
EXP_TYPE == "Bag" ~ "Bag")) %>%
# filter(num_of_euks > 0) %>% # Remove observed zero counts
group_by(SAMPLE, EXPID, EXP_TYPE, IGT_REP, SizeFrac) %>%
summarise(total_euks_noFLP = sum(num_of_euks),
.groups = "rowwise") %>%
# Join with FLP count information
## SAMPLE, EXPID, EXPTYPE, IGTREP, and SizeFrac variables should match
left_join(counts_flp) %>%
replace_na(list(total_FLP = 0, total_euks_wflp = 0)) %>% #Replace NAs with zero
data.frame## Joining, by = c("SAMPLE", "EXPID", "EXP_TYPE", "IGT_REP", "SizeFrac")## Check data frames
# View(counts_flp_compiled)
# View(counts_occur)
# (filter(counts_flp_compiled, SAMPLE == "VD-X18", EXPID == "T40-Rep2"))
# (filter(counts_flp_compiled, SAMPLE == "Piccard-Plume", EXPID == "T10-Rep1"))
# head(counts_flp_compiled)Extract total euk value by adding across nano and micro, then combine back with nano and micro counts.
counts_flp_compiled_all <- counts_flp_compiled %>%
# Exclude size fraction:
group_by(SAMPLE, EXPID, EXP_TYPE, IGT_REP) %>%
summarise(total_euks_noFLP = sum(total_euks_noFLP),
total_FLP = sum(total_FLP),
total_euks_wflp = sum(total_euks_wflp),
.groups = "rowwise") %>%
add_column(SizeFrac = "total_euks") %>% #Add SizeFrac column
bind_rows(counts_flp_compiled) %>% # Combine back with flp compiled list
data.frame
# head(counts_flp_compiled_all)
# filter(counts_flp_compiled_all, SAMPLE == "VD-X18", EXPID == "T40-Rep2")
# filter(counts_flp_compiled_all, SAMPLE == "Piccard-Plume", EXPID == "T10-Rep1")All three size fractions are represented, total, micro, and nano.
# unique(counts_flp_compiled_all$SizeFrac)
# unique(counts_flp_compiled_all$SAMPLE)
# head(counts_flp_compiled_all[1:2,])4.1 Calculate FLP per euk cell calculation
Import metadata which has the extact minutes for each time point used.
metadata <- read.delim("input-data/flp-exp-metadata-compiled.txt")
head(metadata)## SAMPLE TimePoint Minutes SiteOrigin REP
## 1 VD-MustardStand T0 0 Vent Bag
## 2 VD-MustardStand T10 10 Vent Bag
## 3 VD-MustardStand T15 15 Vent Bag
## 4 VD-MustardStand T20 20 Vent Bag
## 5 VD-MustardStand T40 40 Vent Bag
## 6 Piccard-Plume T0 0 Plume BagAdd metadata information to FLP data, reformat sample namse, and calculate FLP per euk as the total FLP divided by the total number of euk cells counted.
counts_flp_calcs_all <- counts_flp_compiled_all %>%
# Add in metadata
# IGTXb are replicate counts, include them as replicates!
separate(EXPID, c("TimePoint", "REP"), sep = "-", remove = FALSE) %>% mutate(
REP = ifelse(grepl("IGT5b", REP), "IGT5", REP),
REP = ifelse(grepl("IGT4b", REP), "IGT4", REP),
REP = ifelse(grepl("Bag", EXP_TYPE), "Bag", REP)) %>%
left_join(metadata, by = c("SAMPLE" = "SAMPLE", "TimePoint" = "TimePoint", "REP" = "REP")) %>%
separate(SAMPLE, c("Site", "Name"), sep = "-", remove = FALSE) %>%
separate(EXPID, c("TimePoint", "Replicate_ID"), sep = "-", remove = FALSE) %>%
## Treat repeated IGT counts completely separate
# group_by(SAMPLE, Site, Name, EXPID, TimePoint, Replicate_ID, EXP_TYPE, IGT_REP, SizeFrac) %>%
## Treat repeated IGT counts as replicates (e.g., IGT4b and IGT4 == IGT4)
group_by(SAMPLE, Site, Name, EXPID, TimePoint, Replicate_ID, EXP_TYPE, REP, SizeFrac) %>%
# FLPperEuk is the total FLP divided by the total number of euk cells counted
mutate(FLPperEuk = total_FLP/(sum(total_euks_noFLP, total_euks_wflp))) %>%
unite("Experiment", Name, REP, sep = "-", remove = FALSE) %>%
data.frame
# filter(counts_flp_calcs_all, SAMPLE == "VD-X18", EXPID == "T40-Rep2")
# filter(counts_flp_calcs_all, SAMPLE == "Piccard-Plume", EXPID == "T10-Rep1")View output
# head(counts_flp_calcs_all[1:2,])
unique(counts_flp_calcs_all$Experiment)## [1] "LotsOShrimp-Bag" "Plume-Bag" "Shrimpocalypse-IGT3"
## [4] "Shrimpocalypse-Bag" "BSW-Bag" "MustardStand-Bag"
## [7] "OMT-IGT4" "Rav2-IGT4" "Rav2-IGT5"
## [10] "Rav2-Bag" "ShrimpHole-Bag" "X18-Bag"# View(counts_flp_calcs_all)COLS: Timepoint, Minutes = time point label, actual incubated minutes COLS: Replicate_ID, REP, and IGT_REP = full replicate identified for IGTs and Bags, designation of biological replicates, and designation of technical replicates for IGT experiments
4.2 Calculate linear regression to obtain slope
Use lm() function in R to calculate linear regression for each experiment. Slope equates to grazing rate. Function inputs the FLP per euk cell data, performs regression and then adds a column for slope and r-squared values.
calculate_lm <- function(df){
regression_1 <- df %>%
type.convert(as.is = TRUE) %>%
## Keep technical replicates separate for IGTs
# group_by(SAMPLE, Site, Experiment, Name, IGT_REP, SizeFrac) %>%
# nest(-SAMPLE, -Site, -Experiment, -Name, -IGT_REP, -SizeFrac) %>%
## Combine technical replicates for IGTs
group_by(SAMPLE, Site, Experiment, Name, REP, SizeFrac) %>%
nest(-SAMPLE, -Site, -Experiment, -Name, -REP, -SizeFrac) %>%
mutate(lm_fit = map(data, ~lm(FLPperEuk ~ Minutes, data = .)),
tidied = map(lm_fit, tidy)) %>%
unnest(tidied) %>%
# select(SAMPLE, Site, Experiment, Name, IGT_REP, SizeFrac, term, estimate) %>%
select(SAMPLE, Site, Experiment, Name, REP, SizeFrac, term, estimate) %>%
pivot_wider(names_from = term, values_from = estimate) %>%
data.frame
# colnames(regression_1) <- c("SAMPLE", "Site", "Experiment", "Name", "IGT_REP",
# "SizeFrac", "INTERCEPT", "SLOPE")
colnames(regression_1) <- c("SAMPLE", "Site",
"Experiment", "Name", "REP",
"SizeFrac", "INTERCEPT", "SLOPE")
out_regression <- df %>%
# group_by(SAMPLE, Site, Experiment, Name, IGT_REP, SizeFrac) %>%
# nest(-SAMPLE, -Site, -Experiment, -Name, -IGT_REP, -SizeFrac) %>%
group_by(SAMPLE, Site, Experiment, Name, REP, SizeFrac) %>%
nest(-SAMPLE, -Site, -Experiment, -Name, -REP, -SizeFrac) %>%
mutate(lm_fit = map(data, ~lm(FLPperEuk ~ Minutes, data = .)),
glanced = map(lm_fit, glance)) %>%
unnest(glanced) %>%
# select(SAMPLE, Site, Experiment, Name, IGT_REP, SizeFrac, r.squared) %>%
select(SAMPLE, Site, Experiment, Name, REP, SizeFrac, r.squared) %>%
right_join(regression_1) %>%
right_join(df) %>%
data.frame
out_regression$SITE <- factor(out_regression$Site, levels = c("VD", "Piccard"))
out_regression$TYPE <- factor(out_regression$EXP_TYPE, levels = c("Bag", "IGT"))
return(out_regression)
}Note that an error may occur when running the below function. This is due to the fact that some experiments did not have replicates.
calcs_wslope_regression <- calculate_lm(counts_flp_calcs_all)## Warning: All elements of `...` must be named.
## Did you want `data = c(EXPID, TimePoint, Replicate_ID, EXP_TYPE, IGT_REP, total_euks_noFLP,
## total_FLP, total_euks_wflp, Minutes, SiteOrigin, FLPperEuk)`?## Warning in summary.lm(x): essentially perfect fit: summary may be unreliable## Warning: All elements of `...` must be named.
## Did you want `data = c(EXPID, TimePoint, Replicate_ID, EXP_TYPE, IGT_REP, total_euks_noFLP,
## total_FLP, total_euks_wflp, Minutes, SiteOrigin, FLPperEuk)`?## Warning in summary.lm(x): essentially perfect fit: summary may be unreliable
## Warning in summary.lm(x): essentially perfect fit: summary may be unreliable## Joining, by = c("SAMPLE", "Site", "Experiment", "Name", "REP", "SizeFrac")
## Joining, by = c("SAMPLE", "Site", "Experiment", "Name", "REP", "SizeFrac")# head(calcs_wslope_regression[1:2,])
# View(calcs_wslope_regression)Use below commented out to recalculate one linear regression. Above function uses the nest() capability of tidyverse. Below, one experiment is subset to checek the value.
# Extract only plume-bag experiment from VD
# tmp_plume <- filter(counts_flp_calcs_all, Experiment == "Plume-Bag") %>% filter(Site == "VD") %>% filter(SizeFrac == "total_euks")
# tmp_plume # View
# Perform linear regression
# lm_out <- lm(FLPperEuk ~ Minutes, data = tmp_plume)
# Check output
# summary(lm_out)
# lm_out$coefficients #Intercept=intercept #Minutes = SLOPE
# Compare with nested function output
# filter(calcs_wslope_regression, Experiment == "Plume-Bag") %>% filter(Site == "VD") %>% filter(SizeFrac == "total_euks") %>% head# View(calcs_wslope_regression)
# head(calcs_wslope_regression)
# unique(calcs_wslope_regression$Site)4.2.1 Plot linear regression trend
calcs_wslope_regression %>%
filter(SizeFrac == "total_euks") %>%
# filter(TYPE != "IGT") %>%
unite(EXPERIMENT, SITE, Experiment, sep = " ", remove = FALSE) %>%
ggplot(aes(x = Minutes, y = FLPperEuk, fill = Site, shape = TYPE, group = Experiment)) +
geom_abline(aes(slope = SLOPE, intercept = INTERCEPT), color = "black", linetype = "dashed", size = 1) +
geom_point(stat = "identity", color = "black",
size = 2, aes(shape = TYPE, fill = Site)) +
scale_shape_manual(values = c(21, 24)) +
scale_fill_manual(values = c("#de2d26", "#1c9099")) +
labs(x = "Minutes", y = bquote("FLP"~eukaryote^-1), title = "Grazing experiment regression") +
facet_wrap(. ~ EXPERIMENT) +
# Report r.squared
geom_text(aes(x = 42, y = max(FLPperEuk), label = paste(round(SLOPE, 4))),
vjust = 1, hjust = 0, size = 3) +
theme_bw() +
theme(strip.background = element_blank(),
strip.text = element_text(color = "black", size = 7),
legend.title = element_blank(),
legend.position = "right")
Data points represent the FLP per euk cells (based on total eukaryote cells counts). Y-axis represents the duration of incubation (in minutes). The dashed purple line reprents the slope and intercept of the experiment.
4.3 Remove Tf-IGT experiments
IGT experiment results appear to have bottle effect, especially in the final time point. Additionally, due to the lack of biological replicates in the IGT experiments, technical replicates are treated as biological replicates in the regression below.
IGT_lm_woTf <- counts_flp_calcs_all %>%
# Select only IGT experiments with total eukaryotes, remove Tf (T3)
filter(SizeFrac == "total_euks") %>%
filter(EXP_TYPE == "IGT" & !(TimePoint == "T3")) %>%
add_column(IGT_cor = "rm Tf") %>%
data.frame
# Recalculate lm(), but keep
igt_regression_noTf <- calculate_lm(IGT_lm_woTf) # Recalculate## Warning: All elements of `...` must be named.
## Did you want `data = c(EXPID, TimePoint, Replicate_ID, EXP_TYPE, IGT_REP, total_euks_noFLP,
## total_FLP, total_euks_wflp, Minutes, SiteOrigin, FLPperEuk,
## IGT_cor)`?
## Warning: All elements of `...` must be named.
## Did you want `data = c(EXPID, TimePoint, Replicate_ID, EXP_TYPE, IGT_REP, total_euks_noFLP,
## total_FLP, total_euks_wflp, Minutes, SiteOrigin, FLPperEuk,
## IGT_cor)`?## Joining, by = c("SAMPLE", "Site", "Experiment", "Name", "REP", "SizeFrac")
## Joining, by = c("SAMPLE", "Site", "Experiment", "Name", "REP", "SizeFrac")Plot IGT grazing experiment with newly calculated grazing effect.
igt_regression_noTf %>%
# filter(SizeFrac == "total_euks") %>%
# filter(TYPE != "IGT") %>%
unite(EXPERIMENT, SITE, Experiment, sep = " ", remove = FALSE) %>%
ggplot(aes(x = Minutes, y = FLPperEuk, fill = Site, shape = TYPE, group = Experiment)) +
geom_abline(aes(slope = SLOPE, intercept = INTERCEPT), color = "black", linetype = "dashed", size = 1) +
geom_point(stat = "identity", color = "black",
size = 2, aes(shape = TYPE, fill = Site)) +
scale_shape_manual(values = c(24)) +
scale_fill_manual(values = c("#de2d26", "#1c9099")) +
labs(x = "Minutes", y = bquote("FLP"~eukaryote^-1), title = "Grazing experiment regression") +
facet_wrap(. ~ EXPERIMENT) +
# Report r.squared
geom_text(aes(x = 5, y = max(FLPperEuk), label = paste(round(SLOPE, 4))),
vjust = 1, hjust = 0, size = 3) +
theme_bw() +
theme(strip.background = element_blank(),
strip.text = element_text(color = "black", size = 7),
legend.title = element_blank(),
legend.position = "right")
4.4 Compile grazing experiment results
calcs_wslope_regression_update <- calcs_wslope_regression %>%
filter(TYPE != "IGT") %>%
bind_rows(igt_regression_noTf %>% select(-IGT_cor)) %>%
data.frame
# Factor
vent_ids <- c("BSW","Plume", "Shrimpocalypse", "LotsOShrimp", "X18", "OMT", "Rav2", "MustardStand", "ShrimpHole")
vent_fullname <- c("Background","Plume", "Shrimpocalypse", "Lots 'O Shrimp", "X-18", "Old Man Tree", "Ravelin #2", "Mustard Stand", "Shrimp Hole")
site_ids <- c("VD", "Piccard")
site_fullname <- c("Von Damm", "Piccard")
# Factor for shipboard
calcs_wslope_regression_update$SiteOrder <- factor(calcs_wslope_regression_update$Site, levels = site_ids, labels = site_fullname)
calcs_wslope_regression_update$NameOrder <- factor(calcs_wslope_regression_update$Name, levels = vent_ids, labels = vent_fullname)
# dim(calcs_wslope_regression); dim(calcs_wslope_regression_update)
head(calcs_wslope_regression_update)## SAMPLE Site Experiment Name REP SizeFrac
## 1 Piccard-LotsOShrimp Piccard LotsOShrimp-Bag LotsOShrimp Bag total_euks
## 2 Piccard-LotsOShrimp Piccard LotsOShrimp-Bag LotsOShrimp Bag total_euks
## 3 Piccard-LotsOShrimp Piccard LotsOShrimp-Bag LotsOShrimp Bag total_euks
## 4 Piccard-LotsOShrimp Piccard LotsOShrimp-Bag LotsOShrimp Bag total_euks
## 5 Piccard-Plume Piccard Plume-Bag Plume Bag total_euks
## 6 Piccard-Plume Piccard Plume-Bag Plume Bag total_euks
## r.squared INTERCEPT SLOPE EXPID TimePoint Replicate_ID EXP_TYPE
## 1 0.66025258 0.3517002 -0.007605829 T0-Rep3 T0 Rep3 Bag
## 2 0.66025258 0.3517002 -0.007605829 T15-Rep3 T15 Rep3 Bag
## 3 0.66025258 0.3517002 -0.007605829 T20-Rep3 T20 Rep3 Bag
## 4 0.66025258 0.3517002 -0.007605829 T40-Rep3 T40 Rep3 Bag
## 5 0.02810998 0.7435088 0.005362856 T0-Rep1 T0 Rep1 Bag
## 6 0.02810998 0.7435088 0.005362856 T0-Rep2 T0 Rep2 Bag
## IGT_REP total_euks_noFLP total_FLP total_euks_wflp Minutes SiteOrigin
## 1 Bag 9 3 2 0 Vent
## 2 Bag 8 4 3 15 Vent
## 3 Bag 9 2 1 20 Vent
## 4 Bag 5 0 0 40 Vent
## 5 Bag 6 6 4 0 Plume
## 6 Bag 5 5 3 0 Plume
## FLPperEuk SITE TYPE SiteOrder NameOrder
## 1 0.2727273 Piccard Bag Piccard Lots 'O Shrimp
## 2 0.3636364 Piccard Bag Piccard Lots 'O Shrimp
## 3 0.2000000 Piccard Bag Piccard Lots 'O Shrimp
## 4 0.0000000 Piccard Bag Piccard Lots 'O Shrimp
## 5 0.6000000 Piccard Bag Piccard Plume
## 6 0.6250000 Piccard Bag Piccard Plume4.5 Check FLP control experiments
All incubations had control experiments run alongside them. This was to ensure added FLP did not decrease or change in concentration over time.
bac_ctrl <- read.delim("input-data/bac-counts-compiled.txt")
# dim(bac_ctrl)
dtaf <- bac_ctrl %>%
separate(SampleID, c("exp", "Replicate", "TimePoint"), sep = "-", remove = FALSE) %>%
separate(Site, c("Site", "Name"), sep = "-", remove = FALSE) %>%
filter(Stain == "DTAF") %>%
data.frame## Warning: Expected 2 pieces. Additional pieces discarded in 17 rows [33, 34, 35,
## 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49].# View(bac_ctrl)
# head(dtaf)
dtaf_avg <- dtaf %>%
group_by(TimePoint, Stain, Site, Name) %>%
summarise(Avg_cellsperml = mean(Cells.ml)) %>%
data.frame## `summarise()` has grouped output by 'TimePoint', 'Stain', 'Site'. You can override using the `.groups` argument.dtaf_avg %>%
filter(Site != "IGT") %>%
ggplot(aes(x = TimePoint, y = Avg_cellsperml, fill = Name, shape = Site)) +
geom_rect(data = filter(dtaf_avg, TimePoint == "T0", Site != "IGT"), aes(
ymin = (Avg_cellsperml-(0.1*Avg_cellsperml)),
ymax = (Avg_cellsperml+(0.1*Avg_cellsperml))), color = NA, alpha = 0.4, xmin = 0, xmax = 6, fill = "black") +
geom_line(aes(group = Name)) +
geom_point(stat = "identity", aes(shape = Site, fill = Name), size = 2) +
# scale_fill_manual(values = c("black","#9970ab", "#5aae61")) +
facet_wrap(Name ~ Site) +
scale_y_log10() +
theme_bw() + theme(strip.background = element_blank(),
legend.title = element_blank(),
axis.text = element_text(size = 10, color = "black"),
title = element_text(size = 10, face = "bold"),
axis.title = element_text(size = 9)) +
labs(title = "FLP counts for controls", y = bquote("Log FLP "~mL^-1), x = "Time point")
dtaf_avg %>%
filter(Site == "IGT") %>%
ggplot(aes(x = TimePoint, y = Avg_cellsperml, fill = Name, shape = Site)) +
geom_rect(data = filter(dtaf_avg, TimePoint == "T0", Site == "IGT"), aes(
ymin = (Avg_cellsperml-(0.1*Avg_cellsperml)),
ymax = (Avg_cellsperml+(0.1*Avg_cellsperml))), color = NA, alpha = 0.4, xmin = 0, xmax = 6, fill = "black") +
geom_line(aes(group = Name)) +
geom_point(stat = "identity", aes(shape = Site, fill = Name), size = 2) +
# scale_fill_manual(values = c("black","#9970ab", "#5aae61")) +
facet_wrap(Name ~ Site) +
scale_y_log10() +
theme_bw() + theme(strip.background = element_blank(),
legend.title = element_blank(),
axis.text = element_text(size = 10, color = "black"),
title = element_text(size = 10, face = "bold"),
axis.title = element_text(size = 9)) +
labs(title = "FLP counts for controls", y = bquote("Log FLP "~mL^-1), x = "Time point")
4.6 Table with grazing values
head(calcs_wslope_regression_update)## SAMPLE Site Experiment Name REP SizeFrac
## 1 Piccard-LotsOShrimp Piccard LotsOShrimp-Bag LotsOShrimp Bag total_euks
## 2 Piccard-LotsOShrimp Piccard LotsOShrimp-Bag LotsOShrimp Bag total_euks
## 3 Piccard-LotsOShrimp Piccard LotsOShrimp-Bag LotsOShrimp Bag total_euks
## 4 Piccard-LotsOShrimp Piccard LotsOShrimp-Bag LotsOShrimp Bag total_euks
## 5 Piccard-Plume Piccard Plume-Bag Plume Bag total_euks
## 6 Piccard-Plume Piccard Plume-Bag Plume Bag total_euks
## r.squared INTERCEPT SLOPE EXPID TimePoint Replicate_ID EXP_TYPE
## 1 0.66025258 0.3517002 -0.007605829 T0-Rep3 T0 Rep3 Bag
## 2 0.66025258 0.3517002 -0.007605829 T15-Rep3 T15 Rep3 Bag
## 3 0.66025258 0.3517002 -0.007605829 T20-Rep3 T20 Rep3 Bag
## 4 0.66025258 0.3517002 -0.007605829 T40-Rep3 T40 Rep3 Bag
## 5 0.02810998 0.7435088 0.005362856 T0-Rep1 T0 Rep1 Bag
## 6 0.02810998 0.7435088 0.005362856 T0-Rep2 T0 Rep2 Bag
## IGT_REP total_euks_noFLP total_FLP total_euks_wflp Minutes SiteOrigin
## 1 Bag 9 3 2 0 Vent
## 2 Bag 8 4 3 15 Vent
## 3 Bag 9 2 1 20 Vent
## 4 Bag 5 0 0 40 Vent
## 5 Bag 6 6 4 0 Plume
## 6 Bag 5 5 3 0 Plume
## FLPperEuk SITE TYPE SiteOrder NameOrder
## 1 0.2727273 Piccard Bag Piccard Lots 'O Shrimp
## 2 0.3636364 Piccard Bag Piccard Lots 'O Shrimp
## 3 0.2000000 Piccard Bag Piccard Lots 'O Shrimp
## 4 0.0000000 Piccard Bag Piccard Lots 'O Shrimp
## 5 0.6000000 Piccard Bag Piccard Plume
## 6 0.6250000 Piccard Bag Piccard Plume# Generate final table
table_grazerate <- calcs_wslope_regression_update %>%
filter(SizeFrac == "total_euks") %>%
group_by(SAMPLE, SiteOrder, NameOrder, SLOPE, EXP_TYPE) %>%
mutate(No_of_Replicates = n_distinct(Replicate_ID)) %>%
select(SAMPLE, Site=SiteOrder, Name=NameOrder, RATE = SLOPE, EXP_TYPE, No_of_Replicates) %>%
distinct() %>%
data.frame
# View(table_grazerate)
# write_delim(table_grazerate, path = "Table-grazingrates.txt", delim = "\t")Add FLP conc to table
# head(table_grazerate)
dtaf_igt <- 5352.8278 # Manually insert FLP concentration for IGT experiments; this value is estimated from how IGT FLP spike-ins were calculated
table_grazerate_wflp <- bac_ctrl %>%
filter(FLP_t0 == "use") %>%
add_column(EXP_TYPE = "Bag") %>%
group_by(Site, EXP_TYPE) %>%
summarise(FLP_conc = mean(Cells.ml)) %>%
right_join(table_grazerate, by = c("Site" = "SAMPLE", "EXP_TYPE" = "EXP_TYPE")) %>%
mutate(FLP_conc = ifelse(EXP_TYPE == "IGT", dtaf_igt, FLP_conc)) %>%
data.frame## `summarise()` has grouped output by 'Site'. You can override using the `.groups` argument.# View(table_grazerate_wflp)Negative grazing rates equal undetected (change to zero).
# head(table_grazerate)
bsw <- c("Plume", "Background")
grazerate_plot <- table_grazerate %>%
type.convert(as.is = TRUE) %>%
mutate(detected = case_when(
RATE < 0 ~ "Not detected",
TRUE ~ "Detected")) %>%
mutate(type = case_when(
Name %in% bsw ~ Name,
TRUE ~ paste("Vent", EXP_TYPE, sep="-")
)) %>%
mutate(GRAZE_RATE = case_when(
RATE < 0 ~ 0,
TRUE ~ RATE
)) %>%
mutate(type_site = case_when(
Name %in% bsw ~ Name,
TRUE ~ "Vent"
)) %>%
# filter(detected == "Detected") %>%
data.frame# Factoring
type_order <- c("Vent-Bag", "Vent-IGT", "Plume", "Background")
grazerate_plot$TYPE <- factor(grazerate_plot$type, levels = type_order)
vent_ids <- c("BSW","Plume", "Shrimpocalypse", "LotsOShrimp", "X18", "OMT", "Rav2", "MustardStand", "ShrimpHole")
vent_fullname <- c("Background","Plume", "Shrimpocalypse", "Lots 'O Shrimp", "X-18", "Old Man Tree", "Ravelin #2", "Mustard Stand", "Shrimp Hole")
site_ids <- c("VD", "Piccard")
site_fullname <- c("Von Damm", "Piccard")
vent_colors <- c("#0868ac", "#41ab5d", "#e7298a", "#c994c7", "#fc4e2a", "#fed976", "#6a51a3", "#ffeda0", "#a1d99b")
names(vent_colors) <- vent_fullname
grazerate_plot$NAME <- factor(grazerate_plot$Name, levels = vent_fullname)# svg("", h =, w = )
grazing_min_plot <- grazerate_plot %>%
ggplot(aes(y = GRAZE_RATE, x = NAME, shape = EXP_TYPE, fill = Site)) +
geom_jitter(stat = "identity", aes(shape = EXP_TYPE, fill = Site),
color = "black", size = 3, width = 0.3) +
scale_shape_manual(values = c(21, 24)) +
scale_fill_manual(values = c("#de2d26", "#1c9099")) +
facet_grid(.~Site, space = "free", scales = "free") +
# coord_flip() +
theme_minimal() +
theme(panel.grid.major = element_line(), panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
axis.text.x = element_text(color="black", size = 12,
angle = 45, hjust = 1, vjust = 1),
axis.text.y = element_text(color="black", size = 12),
axis.title =element_text(color="black", size = 12),
axis.ticks = element_line(),
strip.text =element_blank(), legend.title = element_blank())+
guides(fill = guide_legend(override.aes = list(shape = c(21))),
shape = guide_legend(override.aes = list(fill = "black"))) +
labs(x = "", y = bquote("FLP " ~grazer^-1 ~min^-1))
# dev.off()
grazing_min_plot
5 Plot cell conc & grazing
Generate final plots that show eukaryote and prokaryote biomass and cell concentration. Save final table reporting carbon estimates.
# Subset the average in situ prok cells/ml for non-background samples
tmp <- filter(insitu_proks, Name != "BSW", Name != "Plume") %>% select(MEAN)
avg_insitu <- mean(tmp$MEAN)
# head(insitu_proks)
# Add to master table with data
table_grazerate_wflp_wprok <- insitu_proks %>%
select(Site = SAMPLE, Prok_conc = MEAN, Prok_sem = SEM) %>%
right_join(table_grazerate_wflp) %>%
mutate(Prok_conc = ifelse(is.na(Prok_conc), avg_insitu, Prok_conc)) %>%
data.frame## Joining, by = "Site"Table with all cell count data
# head(table_grazerate_wflp_wprok)
table_grazerate_wflp_wprok_weuk <- plot_euk_format %>%
select(Name = NameOrder, Site.y.y = SiteOrder, euk_conc = avg_conc, EXP_TYPE, euk_conc_sem = SEM_conc) %>%
right_join(table_grazerate_wflp_wprok) %>%
select(SITE = Site.y.y, NAME = Name, EXP = EXP_TYPE, SAMPLE = Site, RATE_min = RATE, REPs = No_of_Replicates, FLP_ml = FLP_conc, PROK_ml = Prok_conc, PROK_sem = Prok_sem, EUK_ml = euk_conc, EUK_sem = euk_conc_sem) %>%
mutate(RATE_min = case_when(
RATE_min < 0 ~ 0,
TRUE ~ RATE_min
)) %>%
data.frame## Joining, by = c("Name", "Site.y.y", "EXP_TYPE")6 Grazing rate calculations
head(table_grazerate_wflp_wprok_weuk)## SITE NAME EXP SAMPLE RATE_min REPs FLP_ml PROK_ml
## 1 Von Damm Background Bag VD-BSW 0.002958889 3 4861.824 37889.62
## 2 Von Damm Plume Bag VD-Plume 0.005274231 3 34242.354 16478.31
## 3 Von Damm X-18 Bag VD-X18 0.001744429 2 3148.722 111429.78
## 4 Von Damm Old Man Tree IGT VD-OMT 0.014510943 2 5352.828 71147.40
## 5 Von Damm Ravelin #2 Bag VD-Rav2 0.003470217 3 51890.942 71147.40
## 6 Von Damm Ravelin #2 IGT VD-Rav2 0.015395240 2 5352.828 71147.40
## PROK_sem EUK_ml EUK_sem
## 1 8608.427 91.83773 21.86613
## 2 2623.935 157.77468 67.09859
## 3 2973.793 314.87222 104.95741
## 4 NA 472.30833 122.45031
## 5 NA 409.33389 73.47019
## 6 NA 620.99799 123.17702Based on Unrein et al. 2007, we use the estimated grazing rate, in situ prok abundance, in situ euk abundance, and the concentration of FLP to make additional estimates.
table_wcalcs <- table_grazerate_wflp_wprok_weuk %>%
# Ingestion rate per hour
mutate(RATE_hr = (RATE_min * 60),
RATE_day = (RATE_hr * 24), #Compare to GR?
# FLP concentration per L
FLP_L = (FLP_ml * 1000),
# mL per grazer per hr
CLEARANCE_RATE_ml = (RATE_hr/FLP_ml),
# nL per grazer per hour
CLEARANCE_RATE_nL = ((RATE_hr/FLP_ml)/1.00E+6),
# proks per grazer per hr
SPEC_GRAZE_RATE_hr = (CLEARANCE_RATE_ml * PROK_ml),
# proks per grazer per day
GRAZE_RATE_DAY = (24 * SPEC_GRAZE_RATE_hr),
# proks per ml per hr
GRAZING_EFFECT_hr = (SPEC_GRAZE_RATE_hr * EUK_ml),
GRAZING_EFFECT_hr_min = (SPEC_GRAZE_RATE_hr * (EUK_ml - EUK_sem)),
GRAZING_EFFECT_hr_max = (SPEC_GRAZE_RATE_hr * (EUK_ml + EUK_sem)),
# cells per ml per day
GRAZING_EFFECT_day = ((SPEC_GRAZE_RATE_hr * 24) * EUK_ml),
# Percentage per day
BAC_TURNOVER_PERC = 100*(GRAZING_EFFECT_day / PROK_ml),
BAC_TURNOVER_PERC_min = 100*(GRAZING_EFFECT_day / (PROK_ml - PROK_sem)),
BAC_TURNOVER_PERC_max = 100*(GRAZING_EFFECT_day / (PROK_ml + PROK_sem))) %>%
data.frame
# View(table_wcalcs)Explaination of units for table with calculated values. * RATE_min & RATE_hr = Grazing rate as ‘FLPs per grazer per minute’ and per hour * CLEARANCE_RATE = ml or nL per grazer per hour * SPEC_GRAZE_RATE (Specific grazing rate) = Prokaryotes per grazer per hour * GRAZING EFFECT = bacteria per ml per hour * Bacterial turnover rate = % per day # Carbon biomass table construction
Unrein, F., Massana, R., Alonso-Sáez, L., and Gasol, J.M. (2007) Significant year-round effect of small mixotrophic flagellates on bacterioplankton in an oligotrophic coastal system. Limnol Oceanogr 52: 456–469.
7 Estimates of carbon biomass and consumption
7.1 Protist carbon-measured
References for estimating biovolume Pernice, M.C., Forn, I., Gomes, A., Lara, E., Alonso-Sáez, L., Arrieta, J.M., et al. (2015) Global abundance of planktonic heterotrophic protists in the deep ocean. ISME J 9: 782–792.
# Import manual biovolume measurements
biov <- read.delim("input-data/biovol-euk-12-10-2020.txt")
# head(biov)
# Calculate volume
biov_calc <- biov %>%
mutate(SizeFrac = case_when(
h >= 20 ~ "micro",
TRUE ~ "nano")) %>%
mutate(Volume = ((pi/6) * (d^2) * d)) %>% # Calculate volume (um cubed) # Hillebrand et al. 1999
mutate(pgC_cell = (0.216 * (Volume^0.939))) %>% # Calculate Cell biomass in pg C per cell # Menden-Deuer and Lessard 2000
data.frame
# View(biov_calc)
biov_calc## EXP VENT_BSW h d SizeFrac Volume pgC_cell
## 1 IGT vent 30.077 25.764 micro 8954.44130 1110.2426245
## 2 IGT vent 89.582 10.000 micro 523.59878 77.1957618
## 3 Bag BSW 14.595 8.036 nano 271.71800 41.6956679
## 4 Bag BSW 12.480 8.982 nano 379.41786 57.0486460
## 5 Bag vent 9.218 3.120 nano 15.90239 2.9015292
## 6 IGT vent 17.255 9.986 nano 521.40274 76.8917043
## 7 Bag vent 41.153 21.000 micro 4849.04826 624.1445904
## 8 IGT vent 10.282 4.136 nano 37.04591 6.4194942
## 9 IGT vent 29.776 25.852 micro 9046.50993 1120.9583343
## 10 IGT vent 10.991 4.000 nano 33.51032 5.8424695
## 11 Bag vent 14.333 2.000 nano 4.18879 0.8290772
## 12 Bag vent 36.164 3.000 micro 14.13717 2.5980292
## 13 Bag BSW 16.206 14.924 nano 1740.42111 238.4669404
## 14 Bag BSW 7.000 7.000 nano 179.59438 28.2640658
## 15 Bag vent 10.069 5.000 nano 65.44985 10.9544849Volume is reported as um^3
# Volume by experiment type
biov_calc %>%
group_by(EXP) %>% summarise(VOL = mean(Volume), C = mean(pgC_cell))## # A tibble: 2 x 3
## EXP VOL C
## * <fct> <dbl> <dbl>
## 1 Bag 836. 112.
## 2 IGT 3186. 400.# Volume by euk size
biov_calc %>%
group_by(SizeFrac) %>% summarise(VOL = mean(Volume), C = mean(pgC_cell))## # A tibble: 2 x 3
## SizeFrac VOL C
## * <chr> <dbl> <dbl>
## 1 micro 4678. 587.
## 2 nano 325. 46.9# Volume by site
biov_calc %>%
group_by(VENT_BSW) %>% summarise(VOL = mean(Volume), C = mean(pgC_cell))## # A tibble: 2 x 3
## VENT_BSW VOL C
## * <fct> <dbl> <dbl>
## 1 BSW 643. 91.4
## 2 vent 2188. 276.# head(biov_calc)
euk_vol <- mean(biov_calc$Volume);euk_vol # in um^3## [1] 1775.759euk_carbon <- mean(biov_calc$pgC_cell); euk_carbon # in pg C per cell## [1] 226.9636euk_carbon_min <- min(biov_calc$pgC_cell); euk_carbon_min## [1] 0.8290772euk_carbon_max <- max(biov_calc$pgC_cell); euk_carbon_max## [1] 1120.958# euk_carbonAvg euk biomass pg C per individual cell == {r}euk_carbon
7.2 Protist carbon-literature
Compare with Menden-Deuer and Lessard 2000, Table 2 - using only the heterotrophic species measured. Based on Table 2, the min volume was 4745 and the maximum was 1.2 x10^7 µm^3. Carbon content was measured at pg per cell, this was 469.48-35,339 pg per cell.
Import the heterotroph species volume and carbon content to compare to my measured values.
# Hu-measured
range(biov_calc$Volume)## [1] 4.18879 9046.50993range(biov_calc$pgC_cell)## [1] 0.8290772 1120.9583343c_prev <- read.delim("input-data/md-lessard-2000.txt") # Table 2, heterotrophs only
# combine and plot
carbon_compare <- c_prev %>%
add_column(source = "Menden-Deuer Lessard") %>%
select(source, Volume = vol, pgC_cell) %>%
rbind(biov_calc %>% add_column(source = "MCR") %>% select(source, Volume, pgC_cell)) %>%
ggplot(aes(x = Volume, y = pgC_cell, fill = source)) +
geom_point(aes(fill = source), shape = 23, color = "black", size = 3) +
scale_y_log10() + scale_x_log10() +
labs(title = "Compare literature to measured cell volume & C content",
x = bquote("Volume" ~µm^-3),
y = bquote("pg C" ~cell^-1)) +
theme_bw() + theme(legend.title = element_blank(),
axis.title = element_text(size = 14),
axis.text = element_text(size = 14),
legend.text = element_text(size = 14))
carbon_compare
Upon comparison, the measured carbon content was much lower from the grazing experiments. This makes sense, as I am looking at preserved specimen and a smaller total number of cells. AND the deep-sea protist cell sizes may be smaller overall.
Find lowest estimates or protist carbon, benthic estimates, and others? How does it compare to my measurements?
euk_carbon_lit_mean <- mean(c_prev$pgC_cell)
euk_carbon_lit_min <- min(c_prev$pgC_cell)
euk_carbon_lit_max <- max(c_prev$pgC_cell)7.3 Prokaryote carbon-literature
Below adding in biomass estimates from prokaryotes and protists.
bac_carbon_ug <- (86)*(1.00E-9) # From Derived from Morono et al. 2011
# bac_carbon_ug
bac_carbon_ug_2 <- (173)*(1.00E-9) # Derived from McNichol et al. 2018; LOFERER-KRO ̈ ßBACHER, J. KLIMA & R. PSENNER 1998
# table_wcalcs7.4 Generate carbon biomass table
Incorporate calculations that include biomass of population and ug C consumed. For rate measurements, only incorporate the Morono et al. 2011 biomass for prokaryotes. This way it is on the lower end and is comparable to Gorda Ridge work.
NEED TO MODIFY SO IT INCLUDES BIOMASS OF EUK AND PROK AT EACH SITE
bsw <- c("Plume", "Background")
table_wcalcs_biomass <- table_wcalcs %>%
add_column(euk_C_ug_Hu = (euk_carbon / (1.00E+06))) %>% # Convert to ug from pg
add_column(euk_C_ug_lit = (euk_carbon_lit_mean / (1.00E+06))) %>% # literature
add_column(bac_C_ug = bac_carbon_ug) %>%
add_column(bac_C_ug_2 = bac_carbon_ug_2) %>%
# Grazing rate in ug C per bac per day
mutate(RATE_ugCbac_pergrazer_perday = (RATE_hr * 24 * bac_C_ug), # Grazing rate as ug C per grazer per day
# % of cell carbon per day
SPEC_INGESTION_RATE = (RATE_ugCbac_pergrazer_perday / euk_C_ug_Hu),
SPEC_INGESTION_RATE_lit = (RATE_ugCbac_pergrazer_perday / euk_C_ug_lit),
Prok_biomass = PROK_ml * bac_carbon_ug,
Euk_biomass_Hu = EUK_ml * euk_C_ug_Hu,
Euk_biomass_lit = EUK_ml * euk_C_ug_lit,
Prok_biomass_L = PROK_ml * bac_carbon_ug * 1000,
Euk_biomass_Hu_L = EUK_ml * euk_C_ug_Hu * 1000,
Euk_biomass_lit_L = EUK_ml * euk_C_ug_lit * 1000,
# Repeat with SEM values
Prok_biomass_sem = PROK_sem * bac_carbon_ug,
Euk_biomass_Hu_sem = EUK_sem * euk_C_ug_Hu,
Euk_biomass_lit_sem = EUK_sem * euk_C_ug_lit,
Prok_biomass_sem_L = PROK_sem * (bac_carbon_ug* 1000),
Euk_biomass_Hu_sem_L = EUK_sem * (euk_C_ug_Hu * 1000),
Euk_biomass_lit_sem_L = EUK_sem * (euk_C_ug_lit * 1000)) %>%
type.convert(as.is = TRUE) %>%
mutate(detected = case_when(
RATE_min < 0 ~ "Not detected",
TRUE ~ "Detected")) %>%
mutate(type = case_when(
NAME %in% bsw ~ NAME,
TRUE ~ paste("Vent", EXP, sep="-")
)) %>%
mutate(GRAZE_RATE = case_when(
RATE_min < 0 ~ 0,
TRUE ~ RATE_min
)) %>%
mutate(type_site = case_when(
NAME %in% bsw ~ NAME,
TRUE ~ "Vent"
)) %>%
data.frame
# View(table_wcalcs_biomass)- Grazing rate column == FLP per minute consumed
- Grazing effect hr == cells per ml per hr
Also make a “bounded” table that demonstrates the ug C consumed in the context of McNichol et al.
# G = number of cells grazed during experiment duration
table_wcalcs_biomass_bounded <- table_wcalcs_biomass %>%
add_column(fgC_cell = 86) %>% # Add in Morono et al. 2011 value
mutate(
# cells_consumed_perday = (G / 1), # Rate of cells consumed * in situ prok, per day
fgC_ml_perday = (GRAZING_EFFECT_day * fgC_cell), # Convert cell amount to fg C
ugC_L_perday = (fgC_ml_perday * (1e-09) * 1000), # Convert to ug C per L
lower_mcnichol = 100*(ugC_L_perday / 17.3),
upper_mcnichol = 100*(ugC_L_perday / 321.4)
) %>%
data.framehead(table_wcalcs_biomass_bounded)## SITE NAME EXP SAMPLE RATE_min REPs FLP_ml PROK_ml
## 1 Von Damm Background Bag VD-BSW 0.002958889 3 4861.824 37889.62
## 2 Von Damm Plume Bag VD-Plume 0.005274231 3 34242.354 16478.31
## 3 Von Damm X-18 Bag VD-X18 0.001744429 2 3148.722 111429.78
## 4 Von Damm Old Man Tree IGT VD-OMT 0.014510943 2 5352.828 71147.40
## 5 Von Damm Ravelin #2 Bag VD-Rav2 0.003470217 3 51890.942 71147.40
## 6 Von Damm Ravelin #2 IGT VD-Rav2 0.015395240 2 5352.828 71147.40
## PROK_sem EUK_ml EUK_sem RATE_hr RATE_day FLP_L CLEARANCE_RATE_ml
## 1 8608.427 91.83773 21.86613 0.1775333 4.260800 4861824 3.651579e-05
## 2 2623.935 157.77468 67.09859 0.3164538 7.594892 34242354 9.241591e-06
## 3 2973.793 314.87222 104.95741 0.1046657 2.511977 3148722 3.324070e-05
## 4 NA 472.30833 122.45031 0.8706566 20.895758 5352828 1.626536e-04
## 5 NA 409.33389 73.47019 0.2082130 4.997113 51890942 4.012512e-06
## 6 NA 620.99799 123.17702 0.9237144 22.169145 5352828 1.725657e-04
## CLEARANCE_RATE_nL SPEC_GRAZE_RATE_hr GRAZE_RATE_DAY GRAZING_EFFECT_hr
## 1 3.651579e-11 1.3835695 33.205668 127.06388
## 2 9.241591e-12 0.1522858 3.654860 24.02685
## 3 3.324070e-11 3.7040034 88.896081 1166.28778
## 4 1.626536e-10 11.5723785 277.737084 5465.73080
## 5 4.012512e-12 0.2854798 6.851515 116.85656
## 6 1.725657e-10 12.2775993 294.662382 7624.36451
## GRAZING_EFFECT_hr_min GRAZING_EFFECT_hr_max GRAZING_EFFECT_day
## 1 96.81058 157.31719 3049.5332
## 2 13.80868 34.24501 576.6444
## 3 777.52519 1555.05037 27990.9067
## 4 4048.68948 6882.77212 131177.5392
## 5 95.88230 137.83081 2804.5574
## 6 6112.04639 9136.68264 182984.7483
## BAC_TURNOVER_PERC BAC_TURNOVER_PERC_min BAC_TURNOVER_PERC_max euk_C_ug_Hu
## 1 8.048465 10.414647 6.558411 0.0002269636
## 2 3.499414 4.162182 3.018725 0.0002269636
## 3 25.119772 25.808540 24.466811 0.0002269636
## 4 184.374334 NA NA 0.0002269636
## 5 3.941897 NA NA 0.0002269636
## 6 257.191065 NA NA 0.0002269636
## euk_C_ug_lit bac_C_ug bac_C_ug_2 RATE_ugCbac_pergrazer_perday
## 1 0.01147298 8.6e-08 1.73e-07 3.664288e-07
## 2 0.01147298 8.6e-08 1.73e-07 6.531607e-07
## 3 0.01147298 8.6e-08 1.73e-07 2.160300e-07
## 4 0.01147298 8.6e-08 1.73e-07 1.797035e-06
## 5 0.01147298 8.6e-08 1.73e-07 4.297517e-07
## 6 0.01147298 8.6e-08 1.73e-07 1.906547e-06
## SPEC_INGESTION_RATE SPEC_INGESTION_RATE_lit Prok_biomass Euk_biomass_Hu
## 1 0.001614483 3.193840e-05 0.003258508 0.02084382
## 2 0.002877822 5.693033e-05 0.001417135 0.03580910
## 3 0.000951827 1.882946e-05 0.009582961 0.07146452
## 4 0.007917726 1.566319e-04 0.006118676 0.10719678
## 5 0.001893483 3.745771e-05 0.006118676 0.09290388
## 6 0.008400232 1.661770e-04 0.006118676 0.14094392
## Euk_biomass_lit Prok_biomass_L Euk_biomass_Hu_L Euk_biomass_lit_L
## 1 1.053653 3.258508 20.84382 1053.653
## 2 1.810146 1.417135 35.80910 1810.146
## 3 3.612524 9.582961 71.46452 3612.524
## 4 5.418786 6.118676 107.19678 5418.786
## 5 4.696281 6.118676 92.90388 4696.281
## 6 7.124700 6.118676 140.94392 7124.700
## Prok_biomass_sem Euk_biomass_Hu_sem Euk_biomass_lit_sem Prok_biomass_sem_L
## 1 0.0007403247 0.004962814 0.2508697 0.7403247
## 2 0.0002256584 0.015228935 0.7698210 0.2256584
## 3 0.0002557462 0.023821507 1.2041746 0.2557462
## 4 NA 0.027791758 1.4048704 NA
## 5 NA 0.016675055 0.8429222 NA
## 6 NA 0.027956696 1.4132080 NA
## Euk_biomass_Hu_sem_L Euk_biomass_lit_sem_L detected type GRAZE_RATE
## 1 4.962814 250.8697 Detected Background 0.002958889
## 2 15.228935 769.8210 Detected Plume 0.005274231
## 3 23.821507 1204.1746 Detected Vent-Bag 0.001744429
## 4 27.791758 1404.8704 Detected Vent-IGT 0.014510943
## 5 16.675055 842.9222 Detected Vent-Bag 0.003470217
## 6 27.956696 1413.2080 Detected Vent-IGT 0.015395240
## type_site fgC_cell fgC_ml_perday ugC_L_perday lower_mcnichol upper_mcnichol
## 1 Background 86 262259.86 0.26225986 1.5159529 0.08159921
## 2 Plume 86 49591.42 0.04959142 0.2866556 0.01542981
## 3 Vent 86 2407217.98 2.40721798 13.9145548 0.74897884
## 4 Vent 86 11281268.37 11.28126837 65.2096438 3.51003994
## 5 Vent 86 241191.93 0.24119193 1.3941730 0.07504416
## 6 Vent 86 15736688.36 15.73668836 90.9635165 4.89629383# View(table_wcalcs_biomass_bounded)
write_delim(table_wcalcs_biomass_bounded, path = "table-wcalc.txt", delim = "\t")
# ?write_delim()7.5 Compilation of plots
Format and factor values to plot. * Grazing rate column == FLP per minute consumed * Grazing effect hr == cells per ml per hr * Specific ingestion rate == % of cell carbon per day, estimated with literature and measured carbon values
biomass_rate_plot <- table_wcalcs_biomass_bounded %>%
select(SITE, NAME, EXP, SAMPLE, type,
PROK_ml, EUK_ml, PROK_sem, EUK_sem,
Prok_biomass_L, Euk_biomass_Hu_L, Euk_biomass_lit_L,
Prok_biomass_sem_L, Euk_biomass_Hu_sem_L, Euk_biomass_lit_sem_L,
GRAZING_EFFECT_hr, GRAZING_EFFECT_hr_min, GRAZING_EFFECT_hr_max, BAC_TURNOVER_PERC, BAC_TURNOVER_PERC_min, BAC_TURNOVER_PERC_max,
GRAZE_RATE, RATE_ugCbac_pergrazer_perday,
SPEC_INGESTION_RATE, SPEC_INGESTION_RATE_lit) %>%
pivot_longer(cols = c(PROK_ml, EUK_ml,
Prok_biomass_L, Euk_biomass_Hu_L, Euk_biomass_lit_L,
GRAZING_EFFECT_hr, BAC_TURNOVER_PERC,
GRAZE_RATE, RATE_ugCbac_pergrazer_perday,
SPEC_INGESTION_RATE, SPEC_INGESTION_RATE_lit),
names_to = "Variable", values_to = "Value") %>%
pivot_longer(cols = c(PROK_sem, EUK_sem, Prok_biomass_sem_L, Euk_biomass_Hu_sem_L, Euk_biomass_lit_sem_L,
BAC_TURNOVER_PERC_max, BAC_TURNOVER_PERC_min, GRAZING_EFFECT_hr_max, GRAZING_EFFECT_hr_min, GRAZING_EFFECT_hr_max),
names_to = "SEM_variable", values_to = "SEM") %>%
data.frame
biomass_rate_plot$NAME_ORDER <- factor(biomass_rate_plot$NAME, levels = c("Background","Plume", "Quakeplume", "Shrimpocalypse", "Lots 'O Shrimp", "X-18", "Old Man Tree", "Ravelin #2", "Mustard Stand", "Shrimp Hole", "Hot Chimlet #1", "Shrimp Gulley", "South of Hot Chimlet", "South of LungSnack", "Arrow Loop", "Bartizan", "Ravelin #1"))
biomass_rate_plot$VARIABLE_ORDER <- factor(biomass_rate_plot$Variable,
levels = c("PROK_ml", "EUK_ml",
"Prok_biomass_L", "Euk_biomass_Hu_L", "Euk_biomass_lit_L",
"GRAZING_EFFECT_hr", "BAC_TURNOVER_PERC",
"GRAZE_RATE", "RATE_ugCbac_pergrazer_perday",
"SPEC_INGESTION_RATE", "SPEC_INGESTION_RATE_lit"),
labels = c("Prokaryote~cells~mL^{-1}", "Eukaryote~cells~mL^{-1}",
"Prokaryote~µg~C~L^{-1}", "Measured~eukaryote~µg~C~L^{-1}", "Literature-based~eukaryote~µg~C~L^{-1}",
"Cells~mL^{-1}~hr^{-1}", "Prokaryote~turnover~'%'~d^{-1}",
"FLP~consumed~min^{-1}", "µg~C~grazer^{-1}~day^{-1}",
"Cell~carbon~%~day^{-1}", "Cell~carbon~%~day^{-1}~(lit)"))
# View(biomass_rate_plot)Function to generate consistent plots for Mid-Cayman Rise cell concentration and biomass data.
conc_rate_plot_mcr <- function(df, var, sem){
df %>%
filter(Variable == var) %>%
filter(SEM_variable == sem) %>%
ggplot(aes(y = Value, x = NAME_ORDER, shape = EXP, fill = SITE)) +
geom_errorbar(aes(ymax = (Value + SEM), ymin = (Value - SEM)),
width = 0.2, position = position_dodge(width = 0.4)) +
geom_point(stat = "identity", aes(shape = EXP, fill = SITE),
color = "black", size = 3, position = position_dodge(width = 0.4)) +
scale_shape_manual(values = c(21, 24)) +
scale_fill_manual(values = c("#de2d26", "#1c9099")) +
facet_wrap(VARIABLE_ORDER ~ ., scales = "free",
strip.position = c("left"), labeller = label_parsed) +
scale_y_log10() +
# scale_y_log10(labels = function(x) format(x, scientific = TRUE)) +
theme_minimal() +
theme(panel.grid.major = element_line(), panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
axis.text.x = element_text(color="black", size = 11,
angle = 45, hjust = 1, vjust = 1),
axis.text.y = element_text(color="black", size = 11),
axis.title =element_text(color="black", size = 14),
axis.ticks = element_line(),
legend.title = element_blank(),
strip.placement = "outside",
strip.text.y = element_text(color="black", size = 11),
strip.text.x = element_blank())+
guides(fill = guide_legend(override.aes = list(shape = c(21))),
shape = guide_legend(override.aes = list(fill = "black"))) +
labs(x = "", y = "")
}# svg("cells-per-ml-yaxies.svg", h = 8, w = 13)
plot_grid(conc_rate_plot_mcr(biomass_rate_plot, "PROK_ml", "PROK_sem"),
conc_rate_plot_mcr(biomass_rate_plot, "EUK_ml", "EUK_sem"),
conc_rate_plot_mcr(biomass_rate_plot, "Prok_biomass_L", "Prok_biomass_sem_L"),
conc_rate_plot_mcr(biomass_rate_plot, "Euk_biomass_Hu_L", "Euk_biomass_Hu_sem_L"))
# dev.off()
# svg("biomass-literature-eukC.svg", h=5, w=7)
# conc_rate_plot_mcr(biomass_rate_plot, "Euk_biomass_lit", "Euk_biomass_lit_sem")
# dev.off()# scale_y_continuous(labels = function(x) format(x, scientific = TRUE)) +
grazing_rate <- biomass_rate_plot %>%
filter(Variable == "GRAZING_EFFECT_hr") %>%
filter(SEM_variable == "GRAZING_EFFECT_hr_min" | SEM_variable == "GRAZING_EFFECT_hr_max") %>%
pivot_wider(names_from = SEM_variable, values_from = SEM) %>%
ggplot(aes(y = Value, x = NAME_ORDER, shape = EXP, fill = SITE)) +
geom_errorbar(aes(ymax = (GRAZING_EFFECT_hr_max), ymin = (GRAZING_EFFECT_hr_min)),
width = 0.2, position = position_dodge(width = 0.4)) +
geom_point(stat = "identity", aes(shape = EXP, fill = SITE),
color = "black", size = 3, position = position_dodge(width = 0.4)) +
scale_shape_manual(values = c(21, 24)) +
scale_fill_manual(values = c("#de2d26", "#1c9099")) +
facet_wrap(VARIABLE_ORDER ~ ., scales = "free",
strip.position = c("left"), labeller = label_parsed) +
theme_minimal() +
theme(panel.grid.major = element_line(), panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
axis.text.x = element_text(color="black", size = 11,
angle = 45, hjust = 1, vjust = 1),
axis.text.y = element_text(color="black", size = 11),
axis.title =element_text(color="black", size = 11),
axis.ticks = element_line(),
legend.title = element_blank(),
strip.placement = "outside",
strip.text.y = element_text(color="black", size = 11),
strip.text.x = element_blank())+
guides(fill = guide_legend(override.aes = list(shape = c(21))),
shape = guide_legend(override.aes = list(fill = "black"))) +
labs(x = "", y = "")
# head(biomass_rate_plot)
# svg("plot-grazing-withscale.svg", h = 7, w = 6)
plot_grid(grazing_rate + scale_y_continuous(limits = c(0,2000)),
grazing_rate,
ncol = 1, rel_heights = c(0.6,1))## Warning: Removed 4 rows containing missing values (geom_point).
# dev.off()
# svg("plot-grazing-log.svg", h = 4, w = 6)
# grazing_rate + scale_y_log10(limits = c(1e-1,1e5), breaks=10^(0:7))
grazing_rate + scale_y_continuous(trans=scales::pseudo_log_trans(base = 10),breaks=10^(0:7))
# dev.off()# scale_y_continuous(labels = function(x) format(x, scientific = TRUE)) +
prok_turnover <- biomass_rate_plot %>%
filter(Variable == "BAC_TURNOVER_PERC") %>%
filter(SEM_variable == "BAC_TURNOVER_PERC_min" | SEM_variable == "BAC_TURNOVER_PERC_max") %>%
pivot_wider(names_from = SEM_variable, values_from = SEM) %>%
ggplot(aes(y = Value, x = NAME_ORDER, fill = SITE)) +
geom_bar(stat = "identity", aes(fill = SITE),
color = "black", width = 2, position = position_dodge2(preserve = "single")) +
geom_errorbar(aes(ymax = (BAC_TURNOVER_PERC_max), ymin = (BAC_TURNOVER_PERC_min)),
width = 0.3, position = position_dodge2(preserve = "single")) +
scale_fill_manual(values = c("#de2d26", "#1c9099")) +
facet_wrap(VARIABLE_ORDER ~ ., scales = "free",
strip.position = c("left"), labeller = label_parsed) +
theme_minimal() +
theme(panel.grid.major = element_line(), panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
axis.text.x = element_text(color="black", size = 11,
angle = 45, hjust = 1, vjust = 1),
axis.text.y = element_text(color="black", size = 11),
axis.title =element_text(color="black", size = 11),
axis.ticks = element_line(),
legend.title = element_blank(),
strip.placement = "outside",
strip.text.y = element_text(color="black", size = 11),
strip.text.x = element_blank())+
guides(fill = guide_legend(override.aes = list(shape = c(21))),
shape = guide_legend(override.aes = list(fill = "black"))) +
labs(x = "", y = "")# svg("prok-turnover.svg", h = 8, w = 6)
plot_grid(prok_turnover,
prok_turnover + scale_y_continuous(limits = c(0,100)),
ncol = 1)## Warning: Removed 3 rows containing missing values (geom_bar).
# dev.off()7.6 Compare MCR data with Gorda Ridge - think about environmental parameters
Import Gorda Ridge data
gr <- read.delim("../../GordaRidgeCruise_2019/protist-gordaridge-2021/Grazing-at-GordaRidge-SKH-2021/data-input/Grazing-calc-wCarbon-results.txt")
# env_gr <- read.delim("../../GordaRidgeCruise_2019/protist-gordaridge-2021/Grazing-at-GordaRidge-SKH-2021/data-input/GR-environ-SAMPLE.txt")
temps <- read.delim("/Users/sarahhu/Desktop/Projects/GordaRidgeCruise_2019/protist-gordaridge-2019-analysis/temperature-allvents.txt")
mcr_graze <- read.delim("table-wcalc.txt")gr## SAMPLE SampleOrigin Bottle Vent.name SAMPLE_ORDER BOTTLE
## 1 BW-Near vent BW BW Exp Near vent BW Near vent BW Experimental
## 2 Vent-Mt Edwards Vent Exp Mt Edwards Mt. Edwards Experimental
## 3 Vent-Venti latte Vent Exp Venti latte Venti latte Experimental
## 4 Vent-Candelabra Vent Exp Candelabra Candelabra Experimental
## 5 Vent-SirVentsalot Vent Exp SirVentsalot Sir Ventsalot Experimental
## cellmL_T0 cellmL_Tf sem_T0 sem_Tf Hrs_Tf Sample.location prok_avg
## 1 69818.95 57808.35 6594.5107 472.550 35 BW 51959.11
## 2 60381.13 40639.10 4786.6367 4095.400 24 Vent 51439.52
## 3 45154.57 32448.30 670.3404 5670.600 19 Vent 111192.50
## 4 25622.60 14701.50 3964.4107 3079.078 26 Vent 55076.66
## 5 27985.33 10606.10 5055.5143 4108.875 18 Vent 52998.29
## MORTALITY MORTALITY_min MORTALITY_max G G_min G_max
## 1 0.1294438 0.06703935 0.1857493 8938.263 4839.402 12329.66
## 2 0.3959460 0.41957541 0.3762200 16818.511 17626.994 16128.80
## 3 0.4174021 0.64113500 0.2325695 31289.007 44259.125 18698.35
## 4 0.5127925 0.57456692 0.4700572 23475.280 25520.903 21977.85
## 5 1.2936683 1.68141744 1.0785052 32912.588 37981.088 29395.13
## GrazingRate_hr GrazingRate_hr_min GrazingRate_hr_max Prok_turnover
## 1 255.3789 138.2686 352.2759 17.20249
## 2 700.7713 734.4581 672.0331 32.69570
## 3 1646.7899 2329.4276 984.1237 28.13949
## 4 902.8954 981.5732 845.3019 42.62292
## 5 1828.4771 2110.0604 1633.0627 62.10122
## Prok_turnover_min Prok_turnover_max N_avg F_avg q G_II_a
## 1 9.313866 23.72954 51959.11 63813.65 0.1882137 9779.414
## 2 34.267415 31.35487 51439.52 50510.12 0.3908531 20105.294
## 3 39.804056 16.81620 111192.50 38801.43 0.3274690 36412.097
## 4 46.337057 39.90411 55076.66 20162.05 0.5416662 29833.162
## 5 71.664736 55.46429 52998.29 19295.72 0.9006783 47734.414
## G_II_b GrazingRate_hr_II fgC_cell_morono fgC_cell_mcnic
## 1 9779.414 279.4118 86 173
## 2 20105.294 837.7206 86 173
## 3 36412.097 1916.4262 86 173
## 4 29833.162 1147.4293 86 173
## 5 47734.414 2651.9119 86 173
## cells_consumed_perday fgC_ml_perday_morono fgC_ml_perday_mcnic
## 1 6129.094 527102.1 1060333
## 2 16818.511 1446391.9 2909602
## 3 39522.957 3398974.3 6837471
## 4 21669.489 1863576.1 3748822
## 5 43883.450 3773976.7 7591837
## ugC_L_perday_morono ugC_L_perday_mcnic lower_mcnichol_morono
## 1 0.5271021 1.060333 3.046833
## 2 1.4463919 2.909602 8.360647
## 3 3.3989743 6.837471 19.647250
## 4 1.8635761 3.748822 10.772116
## 5 3.7739767 7.591837 21.814894
## upper_mcnichol_morono lower_mcnichol_mcnic upper_mcnichol_mcnic
## 1 0.1640019 6.129094 0.3299108
## 2 0.4500286 16.818511 0.9052901
## 3 1.0575527 39.522957 2.1274025
## 4 0.5798308 21.669489 1.1664037
## 5 1.1742305 43.883450 2.3621148all_vents <- mcr_graze %>%
select(SITE, NAME, SAMPLE, EXP, PROK_ml, EUK_ml, GRAZING_EFFECT_hr, GRAZING_EFFECT_hr_min, GRAZING_EFFECT_hr_max, BAC_TURNOVER_PERC, ugC_L_perday) %>%
rbind(gr %>%
add_column(SITE = "Gorda Ridge") %>%
add_column(EUK_ml = NA) %>%
separate(SAMPLE, c("SAMPLE", "NAME"), sep = "-") %>%
select(SITE, NAME, SAMPLE, EXP = Bottle, PROK_ml = prok_avg, EUK_ml, GRAZING_EFFECT_hr = GrazingRate_hr, GRAZING_EFFECT_hr_min = GrazingRate_hr_min, GRAZING_EFFECT_hr_max = GrazingRate_hr_max, BAC_TURNOVER_PERC = Prok_turnover, ugC_L_perday = ugC_L_perday_morono)) %>%
left_join(temps) %>%
mutate(SAMPLE_TYPE = case_when(
grepl("BSW", NAME) ~ "Background",
grepl("Near vent BW", NAME) ~ "Background",
grepl("Background", NAME) ~ "Background",
grepl("Plume", NAME) ~ "Background",
TRUE ~ "Vent"
))## Joining, by = c("NAME", "SAMPLE")# View(all_vents)
all_vents## SITE NAME SAMPLE EXP PROK_ml EUK_ml
## 1 Von Damm Background VD-BSW Bag 37889.62 91.83773
## 2 Von Damm Background VD-BSW Bag 37889.62 91.83773
## 3 Von Damm Background VD-BSW Bag 37889.62 91.83773
## 4 Von Damm Plume VD-Plume Bag 16478.31 157.77468
## 5 Von Damm X-18 VD-X18 Bag 111429.78 314.87222
## 6 Von Damm Old Man Tree VD-OMT IGT 71147.40 472.30833
## 7 Von Damm Ravelin #2 VD-Rav2 Bag 71147.40 409.33389
## 8 Von Damm Ravelin #2 VD-Rav2 IGT 71147.40 620.99799
## 9 Von Damm Ravelin #2 VD-Rav2 IGT 71147.40 620.99799
## 10 Von Damm Mustard Stand VD-MustardStand Bag 56677.00 259.76958
## 11 Von Damm Shrimp Hole VD-ShrimpHole Bag 41982.96 385.71847
## 12 Piccard Plume Piccard-Plume Bag 51429.13 79.30115
## 13 Piccard Shrimpocalypse Piccard-Shrimpocalypse Bag 238585.68 454.81543
## 14 Piccard Shrimpocalypse Piccard-Shrimpocalypse IGT 238585.68 454.81543
## 15 Piccard Lots 'O Shrimp Piccard-LotsOShrimp Bag 53878.14 230.90630
## 16 Piccard Lots 'O Shrimp Piccard-LotsOShrimp Bag 53878.14 230.90630
## 17 Gorda Ridge Near vent BW BW Exp 51959.11 NA
## 18 Gorda Ridge Mt Edwards Vent Exp 51439.52 NA
## 19 Gorda Ridge Venti latte Vent Exp 111192.50 NA
## 20 Gorda Ridge Candelabra Vent Exp 55076.66 NA
## 21 Gorda Ridge SirVentsalot Vent Exp 52998.29 NA
## GRAZING_EFFECT_hr GRAZING_EFFECT_hr_min GRAZING_EFFECT_hr_max
## 1 127.06388 96.81058 157.31719
## 2 127.06388 96.81058 157.31719
## 3 127.06388 96.81058 157.31719
## 4 24.02685 13.80868 34.24501
## 5 1166.28778 777.52519 1555.05037
## 6 5465.73080 4048.68948 6882.77212
## 7 116.85656 95.88230 137.83081
## 8 7624.36451 6112.04639 9136.68264
## 9 793.88963 636.41897 951.36028
## 10 0.00000 0.00000 0.00000
## 11 0.00000 0.00000 0.00000
## 12 44.25930 34.87229 53.64630
## 13 6006.74315 NA NA
## 14 20427.17803 17284.53525 23569.82080
## 15 0.00000 NA NA
## 16 0.00000 NA NA
## 17 255.37893 138.26863 352.27592
## 18 700.77129 734.45809 672.03314
## 19 1646.78986 2329.42764 984.12374
## 20 902.89538 981.57320 845.30189
## 21 1828.47708 2110.06045 1633.06274
## BAC_TURNOVER_PERC ugC_L_perday Highest.Temp pH Mg H2S H2 CH4
## 1 8.048465 0.26225986 4.181 NA NA NA NA NA
## 2 8.048465 0.26225986 4.100 NA NA NA NA NA
## 3 8.048465 0.26225986 4.000 NA NA NA NA NA
## 4 3.499414 0.04959142 4.230 NA NA NA NA NA
## 5 25.119772 2.40721798 62.000 NA NA NA NA NA
## 6 184.374334 11.28126837 121.000 NA NA NA NA NA
## 7 3.941897 0.24119193 112.000 NA NA NA NA NA
## 8 257.191065 15.73668836 112.000 NA NA NA NA NA
## 9 26.780110 1.63858819 112.000 NA NA NA NA NA
## 10 0.000000 0.00000000 108.000 NA NA NA NA NA
## 11 0.000000 0.00000000 22.000 NA NA NA NA NA
## 12 2.065411 0.09135119 4.722 NA NA NA NA NA
## 13 60.423507 12.39791787 90.000 NA NA NA NA NA
## 14 205.482690 42.16169545 90.000 NA NA NA NA NA
## 15 0.000000 0.00000000 36.000 NA NA NA NA NA
## 16 0.000000 0.00000000 35.000 NA NA NA NA NA
## 17 17.202493 0.52710212 NA NA NA NA NA NA
## 18 32.695699 1.44639193 30.000 5.8 42.8 1.01 127.0 10.1
## 19 28.139494 3.39897427 23.000 5.5 50.4 NA NA 0.9
## 20 42.622919 1.86357606 68.000 5.8 45.8 NA 21.9 23.7
## 21 62.101220 3.77397670 53.000 NA 50.8 NA NA NA
## SAMPLE_TYPE
## 1 Background
## 2 Background
## 3 Background
## 4 Background
## 5 Vent
## 6 Vent
## 7 Vent
## 8 Vent
## 9 Vent
## 10 Vent
## 11 Vent
## 12 Background
## 13 Vent
## 14 Vent
## 15 Vent
## 16 Vent
## 17 Background
## 18 Vent
## 19 Vent
## 20 Vent
## 21 VentPlot grazing rates
allrates <- all_vents %>%
select(SITE, NAME, SAMPLE, EXP, SAMPLE_TYPE, starts_with("GRAZING_EFFECT_")) %>%
distinct() %>%
ggplot(aes(y = GRAZING_EFFECT_hr, x = NAME, fill = SITE, shape = SAMPLE_TYPE)) +
geom_errorbar(aes(ymax = (GRAZING_EFFECT_hr_max), ymin = (GRAZING_EFFECT_hr_min)),
width = 0.2, position = position_dodge(width = 0.4)) +
geom_point(stat = "identity", aes(fill = SITE, shape = SAMPLE_TYPE),
color = "black", size = 3, position = position_dodge(width = 0.4)) +
scale_shape_manual(values = c(21, 24)) +
scale_fill_manual(values = c("#de2d26", "#1c9099", "#addd8e")) +
facet_grid(. ~ SAMPLE_TYPE, scales = "free", space = "free") +
theme_minimal() +
theme(panel.grid.major = element_line(), panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
axis.text.x = element_text(color="black", size = 11,
angle = 45, hjust = 1, vjust = 1),
axis.text.y = element_text(color="black", size = 11),
axis.title =element_text(color="black", size = 11),
axis.ticks = element_line(),
legend.title = element_blank(),
strip.placement = "outside",
strip.text.y = element_text(color="black", size = 11),
strip.text.x = element_blank())+
guides(fill = guide_legend(override.aes = list(shape = c(21))),
shape = guide_legend(override.aes = list(fill = "black"))) +
labs(x = "", y = bquote("cells"~ml^-1~hr^-1))# svg("compare-all-rates.svg", h = 10, w = 7)
# plot_grid(allrates,
# allrates + scale_y_continuous(limits = c(0,1000)),
# allrates + scale_y_log10(), ncol = 1)
# dev.off()
# svg("compare-all-rates-color.svg", h = 4, w = 9)
allrates + scale_y_log10()## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Transformation introduced infinite values in continuous y-axis
# dev.off()7.7 Run regression analysis with MCR and GR data
metadata <- read.delim("input-data/flp-exp-metadata-compiled.txt")
# head(metadata)
# temps
all_vents## SITE NAME SAMPLE EXP PROK_ml EUK_ml
## 1 Von Damm Background VD-BSW Bag 37889.62 91.83773
## 2 Von Damm Background VD-BSW Bag 37889.62 91.83773
## 3 Von Damm Background VD-BSW Bag 37889.62 91.83773
## 4 Von Damm Plume VD-Plume Bag 16478.31 157.77468
## 5 Von Damm X-18 VD-X18 Bag 111429.78 314.87222
## 6 Von Damm Old Man Tree VD-OMT IGT 71147.40 472.30833
## 7 Von Damm Ravelin #2 VD-Rav2 Bag 71147.40 409.33389
## 8 Von Damm Ravelin #2 VD-Rav2 IGT 71147.40 620.99799
## 9 Von Damm Ravelin #2 VD-Rav2 IGT 71147.40 620.99799
## 10 Von Damm Mustard Stand VD-MustardStand Bag 56677.00 259.76958
## 11 Von Damm Shrimp Hole VD-ShrimpHole Bag 41982.96 385.71847
## 12 Piccard Plume Piccard-Plume Bag 51429.13 79.30115
## 13 Piccard Shrimpocalypse Piccard-Shrimpocalypse Bag 238585.68 454.81543
## 14 Piccard Shrimpocalypse Piccard-Shrimpocalypse IGT 238585.68 454.81543
## 15 Piccard Lots 'O Shrimp Piccard-LotsOShrimp Bag 53878.14 230.90630
## 16 Piccard Lots 'O Shrimp Piccard-LotsOShrimp Bag 53878.14 230.90630
## 17 Gorda Ridge Near vent BW BW Exp 51959.11 NA
## 18 Gorda Ridge Mt Edwards Vent Exp 51439.52 NA
## 19 Gorda Ridge Venti latte Vent Exp 111192.50 NA
## 20 Gorda Ridge Candelabra Vent Exp 55076.66 NA
## 21 Gorda Ridge SirVentsalot Vent Exp 52998.29 NA
## GRAZING_EFFECT_hr GRAZING_EFFECT_hr_min GRAZING_EFFECT_hr_max
## 1 127.06388 96.81058 157.31719
## 2 127.06388 96.81058 157.31719
## 3 127.06388 96.81058 157.31719
## 4 24.02685 13.80868 34.24501
## 5 1166.28778 777.52519 1555.05037
## 6 5465.73080 4048.68948 6882.77212
## 7 116.85656 95.88230 137.83081
## 8 7624.36451 6112.04639 9136.68264
## 9 793.88963 636.41897 951.36028
## 10 0.00000 0.00000 0.00000
## 11 0.00000 0.00000 0.00000
## 12 44.25930 34.87229 53.64630
## 13 6006.74315 NA NA
## 14 20427.17803 17284.53525 23569.82080
## 15 0.00000 NA NA
## 16 0.00000 NA NA
## 17 255.37893 138.26863 352.27592
## 18 700.77129 734.45809 672.03314
## 19 1646.78986 2329.42764 984.12374
## 20 902.89538 981.57320 845.30189
## 21 1828.47708 2110.06045 1633.06274
## BAC_TURNOVER_PERC ugC_L_perday Highest.Temp pH Mg H2S H2 CH4
## 1 8.048465 0.26225986 4.181 NA NA NA NA NA
## 2 8.048465 0.26225986 4.100 NA NA NA NA NA
## 3 8.048465 0.26225986 4.000 NA NA NA NA NA
## 4 3.499414 0.04959142 4.230 NA NA NA NA NA
## 5 25.119772 2.40721798 62.000 NA NA NA NA NA
## 6 184.374334 11.28126837 121.000 NA NA NA NA NA
## 7 3.941897 0.24119193 112.000 NA NA NA NA NA
## 8 257.191065 15.73668836 112.000 NA NA NA NA NA
## 9 26.780110 1.63858819 112.000 NA NA NA NA NA
## 10 0.000000 0.00000000 108.000 NA NA NA NA NA
## 11 0.000000 0.00000000 22.000 NA NA NA NA NA
## 12 2.065411 0.09135119 4.722 NA NA NA NA NA
## 13 60.423507 12.39791787 90.000 NA NA NA NA NA
## 14 205.482690 42.16169545 90.000 NA NA NA NA NA
## 15 0.000000 0.00000000 36.000 NA NA NA NA NA
## 16 0.000000 0.00000000 35.000 NA NA NA NA NA
## 17 17.202493 0.52710212 NA NA NA NA NA NA
## 18 32.695699 1.44639193 30.000 5.8 42.8 1.01 127.0 10.1
## 19 28.139494 3.39897427 23.000 5.5 50.4 NA NA 0.9
## 20 42.622919 1.86357606 68.000 5.8 45.8 NA 21.9 23.7
## 21 62.101220 3.77397670 53.000 NA 50.8 NA NA NA
## SAMPLE_TYPE
## 1 Background
## 2 Background
## 3 Background
## 4 Background
## 5 Vent
## 6 Vent
## 7 Vent
## 8 Vent
## 9 Vent
## 10 Vent
## 11 Vent
## 12 Background
## 13 Vent
## 14 Vent
## 15 Vent
## 16 Vent
## 17 Background
## 18 Vent
## 19 Vent
## 20 Vent
## 21 VentLinear regression with both Gorda Ridge and MCR data
library(broom)
# ?pivot_longer
regression_input <- all_vents %>%
filter(!is.na(Highest.Temp)) %>%
select(SITE, NAME, SAMPLE, SAMPLE_TYPE, EXP, PROK_ml, EUK_ml, GRAZING_EFFECT_hr, BAC_TURNOVER_PERC, ugC_L_perday, TEMP = Highest.Temp) %>%
pivot_longer(cols = c(GRAZING_EFFECT_hr, BAC_TURNOVER_PERC, ugC_L_perday), names_to = "RATE", values_to = "RATE_VALUE") %>%
pivot_longer(cols = c(PROK_ml, EUK_ml, TEMP), names_to = "PARAMS", values_to = "PARAMS_VALUE") %>%
data.frame
regression_tmp <- regression_input %>%
# Set up the linear regression
group_by(RATE, PARAMS) %>%
nest(-RATE, -PARAMS) %>%
mutate(lm_fit = map(data, ~lm(RATE_VALUE ~ PARAMS_VALUE, data = .)),
tidied = map(lm_fit, tidy)) %>%
unnest(tidied) %>%
select(RATE, PARAMS,
term, estimate) %>%
pivot_wider(names_from = term, values_from = estimate) %>%
select(everything(), SLOPE = PARAMS_VALUE) %>%
data.frame## Warning: All elements of `...` must be named.
## Did you want `data = c(SITE, NAME, SAMPLE, SAMPLE_TYPE, EXP, RATE_VALUE, PARAMS_VALUE)`?regression_results <- regression_input %>%
group_by(RATE, PARAMS) %>%
nest(-RATE, -PARAMS) %>%
mutate(lm_fit = map(data, ~lm(RATE_VALUE ~ PARAMS_VALUE, data = .)),
glanced = map(lm_fit, glance)) %>%
unnest(glanced) %>%
select(RATE, PARAMS, r.squared, adj.r.squared) %>%
right_join(regression_tmp) %>%
right_join(regression_input) %>%
data.frame## Warning: All elements of `...` must be named.
## Did you want `data = c(SITE, NAME, SAMPLE, SAMPLE_TYPE, EXP, RATE_VALUE, PARAMS_VALUE)`?## Joining, by = c("RATE", "PARAMS")
## Joining, by = c("RATE", "PARAMS")# head(regression_results)Plot regression results
regression_results %>%
# filter(RATE == "GRAZING_EFFECT_hr") %>%
# filter(PARAMS == "TEMP") %>%
ggplot(aes(x = PARAMS_VALUE, y = RATE_VALUE, shape = SAMPLE_TYPE, fill = SITE)) +
geom_abline(aes(slope = SLOPE, intercept = `X.Intercept.`), color = "black", linetype = "dashed", size = 0.5) +
geom_point(color = "black", aes(shape = SAMPLE_TYPE, fill = SITE)) +
scale_shape_manual(values = c(21, 24)) +
scale_fill_manual(values = c("#476AA7","#7299CE", "#A2937A")) +
facet_wrap(PARAMS ~ RATE + round(r.squared, 3), scales = "free", ncol = 5,
strip.position = "bottom", labeller = label_parsed) +
theme_bw() +
theme(
strip.background = element_blank(),
strip.placement = "outside",
strip.text = element_text(color = "black", size = 10),
axis.title = element_text(color = "black", size = 10),
legend.title = element_blank()) # +## Warning: Removed 12 rows containing missing values (geom_point).
# labs(y = bquote("Cells "~mL^-1 ~hr^-1), x = "Temperature (C)")sessionInfo()## R version 3.6.1 (2019-07-05)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: macOS Mojave 10.14.6
##
## Matrix products: default
## BLAS/LAPACK: /Users/sarahhu/anaconda3/envs/r_3.6.0/lib/R/lib/libRblas.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] broom_0.7.5 cowplot_1.0.0 forcats_0.5.0 stringr_1.4.0
## [5] dplyr_1.0.4 purrr_0.3.4 readr_1.3.1 tidyr_1.1.3
## [9] tibble_3.1.0 ggplot2_3.3.1 tidyverse_1.3.0
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.1.0 xfun_0.20 haven_2.3.1 colorspace_2.0-0
## [5] vctrs_0.3.6 generics_0.1.0 htmltools_0.5.1.1 yaml_2.2.1
## [9] utf8_1.1.4 blob_1.2.1 rlang_0.4.10 pillar_1.5.0
## [13] withr_2.4.1 glue_1.4.2 DBI_1.1.0 RColorBrewer_1.1-2
## [17] dbplyr_1.4.4 modelr_0.1.8 readxl_1.3.1 lifecycle_1.0.0
## [21] munsell_0.5.0 gtable_0.3.0 cellranger_1.1.0 rvest_0.3.5
## [25] evaluate_0.14 labeling_0.4.2 knitr_1.31 fansi_0.4.2
## [29] highr_0.8 Rcpp_1.0.5 scales_1.1.1 backports_1.2.1
## [33] debugme_1.1.0 jsonlite_1.6.1 farver_2.1.0 fs_1.4.1
## [37] hms_0.5.3 digest_0.6.27 stringi_1.4.6 grid_3.6.1
## [41] cli_2.3.1 tools_3.6.1 magrittr_2.0.1 crayon_1.4.1
## [45] pkgconfig_2.0.3 ellipsis_0.3.1 xml2_1.3.2 reprex_0.3.0
## [49] lubridate_1.7.8 assertthat_0.2.1 rmarkdown_2.6 httr_1.4.2
## [53] rstudioapi_0.11 R6_2.5.0 compiler_3.6.1